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I.  INTRODUCTION 


A.  BACKGROUND 

There  are  many  ways  to  represent  a  three-dimensional  object,  usually  involving 
complicated  mathematical  and  graphical  theories.  Representing  and  displaying  a  three- 
dimensional  object  like  a  terrain  is  important  in  applications  such  as  cartography,  route 
planning  for  robotic  vehicles,  road  construction,  and  planning  of  military  maneuvers. 
A  simplified  model  of  real  terrain  helps  the  user  better  understand  and  analyze  it. 

The  simplest  approximations  are  piecewise  planar,  and  the  very  simplest  is  a 
polyhedron  with  triangular  faces.  Assume  that  a  terrain  is  approximated  by  a  two- 
dimensional  array  of  evenly  spaced  elevations.  Each  group  of  three  adjacent  points  in 
the  array  defines  a  planar  polyhedral  face,  and  a  set  of  triplets  constitutes  a  planar-patch 
approximation  of  the  surface.  Unfortunately,  this  approach  creates  many  small  triangu¬ 
lar  regions,  which  can  make  the  searching  process  of  route  planning  too  complicated. 

In  this  thesis.  I  developed  three  algorithms  for  planar-patch  modeling  of  evenly 
spaced  elevation  data,  and  each  has  its  own  advantages  and  disadvantages.  I  have  used 
the  least-squares  method  to  fit  a  plane  over  the  points  of  a  terrain,  a  quadtree  algorithm 
to  subdivide  regions,  a  gradient  clustering  method  to  build  regions,  and  smoothing  al¬ 
gorithms  on  elevation  data  to  improve  continuity  of  adjoining  regions. 

One  major  advantage  of  planar-patch  modeling  over  other  three-dimensional  surface 
modeling  is  that  terrain  within  a  patch  is  homogeneous  in  the  sense  of  a  constant  gra¬ 
dient.  This  homogeneity  simplifies  calculating  the  cost  or  speed  of  traversing  the  region, 
thus  making  it  easier  to  analyze  for  route  planning. 

Several  thresholds  affect  my  algorithms,  and  different  threshold  values  will  result  in 
different  models.  These  include  the  standard  deviation  of  a  fitted  plane  from  given 
points,  the  ratio  of  magnitude  between  adjacent  points,  and  the  root  mean  square  dif¬ 
ference  of  coefficients  of  planar  equations  between  adjacent  regions.  The  final  result  of 
this  thesis  will  be  provided  to  Ron  Ross  (Ph.D.  student,  Computer  Science  Department, 
Naval  Postgraduate  School)  for  his  research  in  route  planning  for  robotic  vehicles. 

B.  ORGANIZATION 

Chapter  2  introduces  terrain  modeling  techniques,  including  least-squares  planar 
fitting  and  cell  classification. 
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Chapte:  '  is  the  mam  chapter  of  this  thesis,  defining  the  planar-patch  terrain  model. 
Joint  tup-down  and  bottom-up  and  strict  bottom-up  planar-patch  terrain  modeling 
strategies  are  described  in  detail.  Programming  techniques  such  as  the  quadtree  method 
used  in  the  joint  top-down  and  bottom-up  algorithm  and  the  region-growing  method 
used  in  the  strict  bottom-up  algorithm  are  explained.  Finally,  the  issue  of  continuity  of 
adjacent  planar-patches  along  their  common  boundaries  is  discussed.  Attempts  to  solve 
this  problem,  such  as  smoothing  of  original  elevation  data  and  grouping  of  the  same 
type  of  points  into  one  region  are  explained. 

In  Chapter  4  comparisons  are  made  between  the  two  algorithms  in  terms  of  accu¬ 
racy  and  simplicity  of  terrain  representation.  Special  features  and  the  limitations  inher¬ 
ent  to  each  algorithm  are  discussed.  Experimental  results  of  joint  top-down  and 
bottom-up.  strict  bottom-up,  and  presmoothing  bottom-up  algorithms  are  presented  in 
graphical  form  in  this  chapter.  Pictures  of  the  modeled  terrain  are  shown  in  three  di¬ 
mensions  using  pictures  drawn  by  the  Grafstat  and  APL  packages.  Pictures  of  the 
boundaries  of  subregions  are  also  included. 

Finally,  Chapter  5  summarizes  the  contributions  made  by  this  research  and  discusses 
some  of  the  possible  areas  for  further  research  based  on  this  work.  The  Pascal  source 
program  is  in  the  Appendix. 
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II.  REVIEW  OF  TERRAIN  MODELING 


The  most  common  method  of  terrain  representation  is  by  a  set  of  functional  ex¬ 
pressions  such  as 

TWO  =  *3  +  3.vj  -  5y  +>-3 

where  fix.y)  is  the  elevation  function  for  terrain  and 
x  and  y  are  the  latitude  and  longitude  coordinates. 

This  permits  the  simple  calculation  of  the  elevation  given  any  combination  of  x  and  y. 
It  is  known  that  any  continuous  surface  can  b>-  approximated  with  an  arbitrarily  small 
error  by  a  polynomial  of  sufficiently  high  degree.  Such  polynomials  can  be  derived  by 
least-squares  surface-fitting  methods. 

Surfaces  may  be  defined  in  tabular  form  where  the  value  of  f  is  given  for  a  selection 
of  representative  arguments.  To  find  an  arbitrary  value  of  f,  a  table  lookup  can  be  per¬ 
formed  to  determine  the  value  of  the  nearest  known  points  and  use  them  to  approximate 
the  desired  value  by  interpolation.  [Ref.  1:  p.  212] 

It  is  possible  to  specify  a  surface  in  terms  of  guiding  points  by  a  generalization  of 
the  Bezier  polynomials  or  the  B-splines.  The  Bezier  polynomial  or  B-splinc  method  finds 
an  approximate  curve  that  passes  near  a  set  of  given  points.  [Ref.  2:  p.  247] 

A.  TERRAIN  CELL  CLASSIFICATION 

1.  QUADRATIC  APPROXIMATION  OF  SURFACE  BY  LEAST  SQUARES 

FIT 

Fitting  a  least-squares  quadratic  to  a  3  by  3  square  of  points  on  a  surface  is  done  in  the 
following  way.  The  fit  as  a  function  of  x  and  y  is:  [Ref.  3:  p.  55] 

fixyj  =  kl  +  k2x  +  kxv  +  kAx  +  k5xy  +  k^y 

The  square  of  the  error  between  the  fit  equation  and  the  given  terrain  points  is: 

E1  =  +  k*x  +  k^'  +  *4 x'ksxy  +  k^y2  -fxy))2 

*  y 

The  partial  derivatives  of  the  squared  error  with  respect  to  each  constant  are: 


Yj-  =  +  An*  +  kxy  4-  A4.v2A5.rv  +  Aav2  —  ./Uv))  1 

*  L  .X  V 


£=[£!> 

1  L  x  y 


2  2 

i  +  k2.x  +  kxv  +  A4. x  ksxy  +  kxv  -f[xy))  x 


~ ~j\ —  =  2^* '  V  /^i  "b  Avv  +  A3V  +  A4.v  k$xy  +  k^y' 
3  L  jc  >■ 


-ftxy))  y 


#4l2> 


i  +  kjX  +  kxv  -f  k4x2ksxy  +  k^y2 


-7(-vj'))  ■> 


=  2£5>.  +  k2x  +  A'3_v  +  k4x2ksxy  +  A^v2  -J[xy))  x 

5  L  A  ‘  V J 


~  =  2 V V(A,  4-  A2*  +  kxv  +  k4x2ksxy  +  k<y2 -J{xy))  y2  (10) 

6  La  v 

To  minimize  the  squared  error,  each  of  the  above  equations  is  set  to  zero  and  simplified. 
Those  equations  can  be  further  simplified  by  assuming  that  the  coordinates  of  each  of 
the  points  in  the  3  by  3  matrix,  except  the  center  point,  is  ±  one  distance  unit  from  the 
center  point.  By  the  above  assumption,  any  terms  in  above  equations  which  contain  a 
summation  of  an  x  or  y  coordinate  with  an  odd  power  will  go  to  zero.  Simplified 
equations  are: 


l=A-,yy  +  A4x6  +  A6x6  —  22m 


0  =  A,  x  6  4- 


0  =  A,  x  6  + 


22m tv 


4 


(14) 


0  =  A,  x  6  +  A4  x  6  +  A6  x  4  -  y  Yjlx^x2 

x  y 

0  =  A;X4  -  Y^AX'}'  )'r-V'  (  1  5) 

x  y 

0  =  A,  x  5  +  A4  x  4  +  k6  x  6  -  (16) 

x  y 

From  above  equations,  if  we  solve  for  A,.  k2,  k},  A„.  ks  and  A6,  we  get  a  quadratic 
equation  fitted  to  the  3  by  3  square  of  points. 

2.  TERRAIN  CELL  CLASSIFICATION  BASED  ON  QUADRATIC 
APPROXIMATIONS 

The  eigenvalues  of  the  Hessian  matrix  for  the  quadratic  function  approximated 
for  the  3  by  3  points  on  the  surface  are:  [Ref.  3:  p.  61) 


(2k.  +  2k6)  ±  N  (  -2k,  -  2 k0Y  -  4(2A42A9  -  kj) 


Using  the  signs  of  the  eigenvalues,  we  can  classify  the  central  point  by  Table  1 
on  page  6.  [Ref.  3:  p.  69] 
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Table  1.  CLASSIFICATION  OF  SURFACE 
POINTS  BY  EIGENVALUES  OF  THE 
HESSIAN  MATRIX 


sign  of 

sign  of  ).2 

- 

0 

+ 

+ 

saddle 

impossible 

depression 

0 

ridge 

plane 

valley 

- 

hill 

impossible 

pass 

B.  PLANE  APPROXIMATION  OF  A  SURFACE  BY  LEAST  SQUARES  FIT 

The  elevation  z  of  any  point  (x.y)  can  he  expressed  as: 

r  =A*,v) 


The  plane  that  we  want  to  fit  over  the  surface  is: 


z  =  ax  -r  by  4-  c 

The  square  of  the  error  between  the  fitted  plane  and  the  real  terrain  is: 
E1  =  XE<«  +  by  +  c  -J[xy))2 

X  >• 

The  partial  derivatives  of  the  squared  error  with  respect  to  a,b  and  c  are: 

•-£-  =  [2ZZ(ajc  +  by+c -A**))  ]* 


(21) 
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cEl 

cb 


2Vy^(ar  4  by  4  c  -J[xy))  [>■ 

X  V  J 


(22) 


=  2^^(ov  4  by  4  c  -j[xy))  1 

CC  L  x  y 


(23) 


To  minimize  the  squared  error,  each  of  the  above  equations  is  expanded  and  then  the 
partials  are  set  to  zero. 


° = xiX/u'2 + YjYjbxy + 1lEjcx  ~ 


(2d) 


xv  x  y 


X  v  X  y 


o  -  ZZ">- + ZZ*: + ZZc>-  - 


(25) 


A  V 


YjYjhy + XX/  ~  XXX*^ 


(26) 


Here,  we  have  three  unknowns,  a,b  and  c,  and  three  equations.  So  we  can  solve  for 
them. 

C.  THE  GRADIENT  AT  A  POINT 

One  of  the  primary  criterion  we  will  use  for  grouping  points  into  a  region  is  the 
gradient.  The  gradient  of  a  point  is  a  vector  that  has  a  magnitude  (slope)  and  a  direction 
for  its  components.  We  can  represent  a  three  by  three  elevation  matrix  of  terrain  points 
as  in  Table  2  on  page  8,  where  x  and  y  are  the  coordinates  of  a  point  on  a  terrain  and 
f^x,y)  is  the  elevation  of  the  point. 
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Table  2.  A  MATRIX  OF  TERRAIN 


POINTS 


f(x-l,y+  1) 

Hx,y+I) 

fl[x+l,v+l) 

f(x-l,y) 

fU,y) 

Hx+I,y) 

fTx-l.y-1) 

«x,y-I) 

Rx+  l.y-l) 

The  gradient  vector  of  a  point  (x.y)  can  be  approximated  in  three  ways.  The  first 
method  is: 

Ai  =AX  +  lj) 

a2  =Av+  i) -A*S) 

magniiude{. rvv)  =  ^  (A;  +  At) 


-1  a2 

direciion(xy)  =  tan  ( — — ) 


The  second  method  is: 

Ai  =A*+  l>y)-Ax-  l»v) 
a2  =/0u'  +  i) -Ax<y-  i) 

The  magnitude(.x,y)  and  direction(x.y)  are  the  same  as  above. 
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The  third  way  of  approximating  gradient  is  by  first  fitting  a  quadratic  to  the  three 
by  three  square  of  points  and  then  using  the  coefficients  of  the  quadratic  :  [Ref.  3:  p.  60] 


magnitude^. rj)  =  N  (Af  +  k2) 


direction(xy)  —  tan  ( -j—  ) 
*2 


where  k2  and  /c3  are  the  coefficients  of  the  quadratic  described  in  section  A. 

The  first  method  considers  only  one  quadrant  of  the  3  by  3  matrix,  thus  is  biased. 
The  second  method  uses  a  point  from  each  of  the  quadrants,  and  is  better  than  the  first 
method  in  the  sense  that  it  utilizes  more  unbiased  information.  However,  it  gives  an 
incorrect  value  in  the  case  as  shown  in  Table  3. 


Table  3.  AN  EXAMPLE  OF  A  TER¬ 
RAIN  POINTS  MATRIX 


5 

2 

3 

1 

3 

2 

4 

When  we  calculate  the  gradient  of  center  point  using  the  second  method,  it  will  give 
the  magnitude  value  zero  and  direction  value  45  degree,  which  is  wrong.  The  third 
method  does  not  have  such  problem,  since  it  utilizes  all  eight  neighboring  points.  In  this 
case  it  gives  the  magnitude  value  0.2357  and  direction  value  45  degree.  I  used  the  third 
method  in  my  program. 
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The  gradient  vector  of  a  terrain  point  points  in  the  x-y  direction  of  maximum  slope 
of  that  point;  magnitude' xy)  is  the  maximum  slope  value  and  direction f xjy  is  the  di¬ 
rection  of  that  slope  in  the  x-y  plane. 
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III.  PLANAR-PATCH  TERRAIN  MODELS 


A  planar-patch  terrain  model  is  a  model  that  represents  terrain  in  the  form  of  a 
polyhedron.  This  model  is  different  from  the  common  planar  terrain  model  used  in 
graphics  in  that  size  of  a  planar-patch  is  not  fixed,  and  the  shape  of  a  patch  does  not 
have  to  be  regular. 

A.  IMPLEMENTATION 

I  chose  Pascal  as  an  implementation  language  for  this  research  because  of  my  fa¬ 
miliarity  and  experience,  and  the  Waterloo  Pascal  on  the  IBM  370  3033AP  mainframe 
was  the  Pascal  compiler  that  was  used. 

The  APL  language  and  the  Grafstat  graphics  software  package  were  used  to  verify 
the  correctness  and  accuracy  of  the  algorithms.  Even  though  I  chose  Pascal  as  my  im¬ 
plementation  language,  the  matrix  form  of  the  elevation  data  suggests  APL.  With  the 
Grafstat  software  package  which  runs  in  the  APL  environment,  one  can  check  the 
intermediate  results  of  an  APL  program  by  displaying  pictures.  Since  I  experimented 
with  different  threshold  values  in  my  algorithm.  APL  and  Grafstat  were  useful.  The 
three-dimensional  pictures  in  chapter  4  were  done  with  Grafstat. 

Two  arrays  of  records  are  used  as  main  data  structures  by  all  my  algorithms: 
point_array  and  subregion_array.  The  point_array  is  a  two-dimensional  array  which 
corresponds  to  the  points  in  the  real  terrain.  Each  element  of  the  array  is  a  record  which 
has  the  elevation,  the  subregion  ID.  and  the  magnitude  of  the  terrain  point  as  its  fields. 
From  the  values  of  elevation  and  magnitude,  the  subregion  ID  is  determined  and  as¬ 
signed  to  the  point.  The  subregion_array  is  a  one-dimensional  array  of  records  which  is 
created  during  the  subregion  creation  phase.  The  indices  of  this  array  are  the  ID  s  of 
subregions.  It  keeps  necessary  information  about  the  subregions  such  as  the  equation 
of  the  plane  fitted  to  the  subregion,  the  boundary  of  the  subregion,  and  other  subsidiary7 
information. 

The  general  structure  of  my  first  two  modeling  algorithms  is  shown  in  Figure  1  on 
page  12. 
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B.  JOINT  TOP-DOWN  AND  BOTTOM-UP  PLANAR-PATCH  TERRAIN 
MODELING 

My  first  algorithm  consists  of  two  parts.  The  Hrst  part  ic  the  top-down  subdivision 
of  the  original  elevation  data  matrix  into  subregions  and  the  second  part  is  the 
bottom-up  merging  of  adjoining  subregions  when  adjoining  subregions  have  similar 
plane  equations. 
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1.  TOP-DOWN  PHASE 


A  quadtree  region  subdivision  method  is  used  during  the  top-down  phase  of  the 
firtst  algorithm.  The  method  generates  a  quadtree  by  recursively  dividing  a  two- 
dimensional  region  into  quadrants.  Each  node  in  the  quadtree  has  four  data  elements, 
one  for  each  of  the  quadrants  in  the  region.  [Ref.  4:  p.  215]  The  original  region  is  given 
some  standard  tests  such  as  comparison  to  a  threshold  of  the  maximum  magnitude  or 
the  minimum  magnitude  of  the  points  within  it,  depending  on  the  application.  Only  if 
a  region  fails  the  tests  and  is  not  of  minimum  size,  it  is  further  subdivided  into  four  su¬ 
bregions.  If  a  region  has  both  length  and  width  four  cells  or  more,  it  is  not  of  the  min¬ 
imum  size.  See  Figure  4  on  page  16  for  the  different  kinds  of  minimal  subregions.  The 
subdivision  process  described  above  is  repeated  for  each  subregion.  Figure  2  on  page 
14  shows  graphically  this  algorithm. 

The  code  'xorks  whether  a  region  has  a  dimensions  of  odd  or  even  numbers.  For 
example,  consider  4  by  4  and  5  by  5  matrix  regions  to  subdivide;  Figure  4  on  page  16 
shows  how  this  method  works. 

2.  THRESHOLDS  USED 

Three  thresholds  are  used  in  the  top-down  phase:  low_bound,  upper_bound  and 
standard_deviation  (SD). 

The  threshold  upper_bound  is  set  by  the  maximum  climbing  capability  of  the 
vehicle  in  cases  where  the  polyhedral  model  will  be  used  for  route  planning.  If  the  ab¬ 
solute  value  of  the  minimum  slope  of  a  region  is  greater  than  the  upper_bound.  which 
means  every  cell  has  a  slope  greater  than  the  upper_bound.  th^n  the  region  is  assumed 
being  untraversable  by  the  vehicle,  and  the  region  need  not  be  subdivided  further  be¬ 
cause  a  good  polyhedral  fit  just  doesn  t  matter  for  that  region.  The  threshold  low_bound 
means  that  when  the  absolute  value  of  the  maximum  slope  of  a  region  is  less  than  the 
low_bound.  which  means  every  cell  has  a  slope  less  than  the  low_bound.  then  the  slopes 
in  the  region  are  unimportant.  Such  regions  need  not  be  subdivided  either. 

For  the  region  which  has  maximum  slope  greater  than  the  low_bound  and  mini¬ 
mum  slope  less  than  the  upper_bound,  the  region  is  subdivided  in  four  parts  according 
to  the  quadtree  algorithm  to  better  isolate  the  unlevel  slopes.  Then  we  fit  a  plane  to  the 
region  and  calculate  the  standard  deviation  of  the  plane  fom  the  real  terrain  points. 

The  standard  deviation  (SD)  is: 
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PROCEDURE  MAKE_SUBREGIONS( YL , YU. XL , XU  :  INTEGER; 

prev_region  :  region_p,tr) ; 

(X  YL  :  LOWER  BOUND  IN  Y  AXIS,  YU  :  UPPER  BOUND  IN  Y  AXIS, 

XL  :  LOWER  BOUND  IN  X  AXIS,  XU  ••  UPPER  BOUND  IN  X  AXIS  X) 

VAR 

P  :  REGION.PTR; 

BEGIN 

IF  ((YU-YL)  >  2)  AND  ((XU-XL)  >  2)  THEN  (X  NOT  MINIMUM  SIZE  *) 

if  C p3 . max_mag_of_region  >  low_bound)  then  (X  test  x) 
begin 

make_subregions( ( (yu-yl )  div  2)+yl+l  ,  yu, 

((xu-xl)  div  2)+xl+l  ,  xu  ,  p);  (x  1st  quadrant  x) 
make_subregions(yl  ,  ((yu-yl)  div  2)+yl, 

((xu-xl)  div  2)+xl+l  ,  xu  ,  p);  (x  2nd  quadrant  X) 
make_subregions(yl  ,  ((yu-yl)  div  2)+yl, 

xl  ,  ((xu-xl)  div  2)+xl,  p);  (*  3rd  quadrant  K) 

make_subregions( ( (yu-yl )  div  2)+yl+l  ,  yu, 

xl  ,  ((xu-xl)  div  2)+xl  ,  p);  (*  4th  quadrant  *) 

P  :=  prev_region  ;  (*  goto  parent  region  *) 

end 

else  stop_subdi vide  («  region  passed  standard  test  K) 
else  stop_subdi vide  (x  region  has  minimum  size  *) 

END 


Figure  3.  Pascal  Code  for  the  Quadtree  Region  Subdivision  Method 


where  J[xj)  is  the  elevation  of  point  (x,y) 

ax  +  by  *F  c  is  the  plane  equation  of  the  region 
n  is  the  total  number  of  points  in  that  region. 

The  region  is  subdivided  only  if  the  calculated  SD  is  bigger  than  the  threshold 
SD  and  it  does  not  have  the  minimum  size. 

3.  BOTTOM-UP  PHASE 

We  now  merge  adjoining  regions  having  similar  plane  expressions.  The  reason 
is  that  in  the  process  of  applying  the  quadtree  subdivision,  we  often  subdivide  a  region 
which  should  not  have  been  so  completely  subdivided.  For  example,  say  a  region  has 
very  steep  slopes  in  the  first  quadrant  and  the  rest  of  the  quadrants  are  flat.  Since  the 
original  region  will  fail  in  the  maximum  slope  test  due  to  the  slopes  in  the  first  quadrant, 


15 


Figure  A.  Minimum-Size  Subregions  Created  by  the  Quadtree  Method:  This  pic 
ture  shows  quadtree  subdivision  method  applied  to  the  regions  which 
have  dimensions  of  odd  or  even  numbers,  and  also  shows  boundaries 
and  different  kinds  of  minimum  si/e  subregions. 

the  top-down  phase  will  subdivide  it  in  four.  However,  the  first,  second  and  third 
quadrants  are  all  flat  and  should  not  be  separate. 

The  adjoining  regions  of  a  region  arc  found  from  the  point_array  by  looking  up 
the  subregion  ID  field  of  adjoining  point.  The  bottom-up  phase  is  done  in  two  sub- 
phases:  row_contbine  and  col_combine.  The  ro\»_cotubine  finds  and  merges  adjoining  re¬ 
gions  which  are  located  in  a  same  row,  and  the  col_comnine  finds  and  merges  the  regions 
in  a  same  column. 
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Figure  5.  Subregions  after  the  Merging  Process:  These  subregions  are  created 
by  performing  the  merging  process  o\cr  the  initial  subregions  created  by 
the  quadtree  algorithm. 

To  merge  adjoining  subregions  which  have  similar  plane  equations,  a  value 
equation_difference  is  calculated  from  the  two  plane  equations  and  compared  to  the 
threshold  abc_difference.  The  equation_differenee  is: 

equation  jliffercnce  =  ^  (a.  —  a2)‘  -r  (b]  —  b2)  +  (Cj  —  c2)‘ 

where  at,  bt  and  c,  are  the  coefficients  ot  plane  x  and 
a,.  b3  and  c,  are  the  coefficients  of  plane2- 
If  this  is  less  than  the  threshold  then  the  two  adjoining  regions  pass  the  test. 
Then  we  fit  a  plane  over  the  two  subregions  and  calculate  the  standard  deviation  of  the 


plane,  and  only  if  it  is  less  than  the  threshold  standard_deviation  we  do  merge  the  two 
subregions. 

C.  STRICT  BOTTOM-LP  PLANAR-PATCH  TERRAIN  MODELING 

1.  POINT-CLASSIFICATION  PHASE 

In  strict  bottom-up  modeling,  we  first  classify  every  point  in  a  region  by  at¬ 
taching  two  tags  to  it.  The  first  tag  classifies  slope  of  the  point:  level,  safe-slope  or 
unsafe-slope.  A  threshold  lo\v_bound  defines  the  low  limit  of  slope  and  another  threshold 
upper_b°und  defines  the  upper  limit.  The  second  tag  classifies  the  curvature  of  the  point, 
by  the  terrain  cell  classification  criteria  described  in  chapter  2.  The  designators  are  hill, 
depression,  plane  and  special.  The  tag  special  refers  to  all  the  rest  of  the  cell  classifica¬ 
tions  except  the  above  three.  A  threshold  called  curvature  is  compared  to  the  eigenvalue 
of  the  Hessian  matrix;  if  its  absolute  value  is  less  than  the  threshold,  then  the  eigenvalue 
is  regarded  as  zero. [Ref.  3:  p.  69]  We  group  the  contuguous  points  with  the  same  first 
and  second  tag  to  form  subregions. 

Next  the  thresholds  magnitude_ratio  and  abc_difference  are  used  to  combine  the 
subregions.  The  magnitude_ratio,  is  used  in  the  region-growing  phase  and  the 
abe_difference  is  used  as  in  joint  top-down  and  mottom-up  modeling. 

2.  REGION-GROWING  PHASE 

The  region-growing  phase  consists  of  two  subphases:  one-dimensional  and 
two-dimensional. 

«.  ONE-DIA I ENSIONA L  REGION- GROW 7  AG 

'Ihe  one-dimensional  phase  starts  with  the  first  point  (pn)  in  the  first  row 
of  the  point_array.  The  point  becomes  the  head  of  a  linked  list  and  the  subregion_ID 
field  of  the  point  is  given  the  value  one.  The  subregion  ID  starts  from  one  and  is  in¬ 
cremented  every  time  a  new  subregion  is  created.  Whenever  a  new  subregion  is  created, 
the  necessary  information  about  the  subregion,  such  as  the  coordinates  of  the  starting 
point  and  boolean  value  that  says  whether  the  subregion  is  active  or  not,  is  stored  in  the 
subregion_array  using  the  subregion  ID  as  an  index.  Then  it  examines  the  next  point  ( 

p,,)  in  the  same  row.  If  the  ratio  computed  using  magnitudes  of  gradient  ( 
abs{magnitudepV,  —  magnitude ,,) 

- — - — - — - )  is  less  than  the  threshold  magnitude  ratio,  then  they 

tnci^Yiittidc | 

are  linked  together  to  form  a  linear  subregion  (Ru).  That  is,  the  next_point  field  of  the 
first  point  in  the  point_array  is  given  the  x  and  y  coordinate  of  the  second  point  and  the 
subregion_ID  field  of  the  second  point  is  given  the  same  value  as  the  subregion  ID  of  the 
first  point.  The  average  magnitude  of  the  gradient  of  the  points  within  the  subregion  is 
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Figure  6.  One-Dimensional  Region-Growing 


calculated  and  is  stored  in  the  avg_mag  Field  of  the  subregion  record  in  the 
subregion_array.  Then  repeatedly  for  each  next  point  (/>„')  in  the  same  row.  the  ratio 
between  the  average  magnitude  of  tiie  gradient  of  the  subregion  and  the  magnitude  of 
tire  gradient  of  the  new  point  is  compared  to  magnitude_ratio.  If  on  this  basis  new  point 
(;>..)  should  not  join  the  subregion,  a  new  region  starts  from  that  new  point  and  the  new 
point  becomes  a  head  of  another  linked  list.  This  process  continues  for  all  rows.  See 
Figure  6  for  an  example  of  this  algorithm. 

b.  TWO-DIMENSIONAL  REG  I OK-G  ROWING 

The  two-dimensional  subphase  starts  with  the  first  linear  subregion  (/?,,) 
which  was  created  by  the  one-dimensional  subphase.  Actually,  it  starts  from  the  first 
point  in  the  first  column  of  the  point_array.  which  belongs  to  the  subregion  one.  It  ex¬ 
amines  the  point  below  in  the  same  column,  which  belongs  to  another  subregion,  to 
decide  whether  the  average  magnitude  of  the  gradient  of  the  subregions  are  similar.  The 
information,  average  magnitudes  of  the  gradient  of  the  two  subregions,  is  retrieved  from 
the  the  subregion_array  using  the  subregion  ID  of  each  point.  If  the  ratio  is  less  than 
magnitude_ratio  then  they  are  linked  together.  The  tail  of  the  first  subregion  is  linked 


19 


Figure  7.  Two-Dimensional  Region-Growing:  This  algorithm  is  applied  after  the 
one-dimensional  region-growing  algorithm. 


to  the  head  of  the  second  subregion,  the  tail  of  the  second  subregion  is  left  nil.  all  the 
subregion  IDs  of  the  points  within  the  second  subregion  are  changed  to  the  first  subre¬ 
gion  ID,  and  the  second  subregion  in  the  subregion_array  is  marked  inactive.  When 
those  two  linear  subregions  arc  merged,  the  resulting  subregion  is  not  one-dimensional 
any  more,  except  if  both  subregions  were  single  points.  Then  repeatedly  for  each  next 
point  in  the  same  column,  the  ratio  between  the  average  magnitudes  of  the  gradient  of 
the  subregions  is  compared  to  the  magnitude_ratio.  If  a  new  subregion  does  not  join  the 
old  subregion,  the  control  is  moved  to  the  new'  subregion  and  the  process  starts  from 
that  new  subregion.  This  process  continues  for  all  columns.  This  is  illustrated  in 
Figure  7.  After  the  two-dimensional  phase,  we  fit  a  plane  to  each  subregion  using 
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least-squares.  Then  the  threshold  abc_difference  is  used  to  merge  adjoining  regions  that 
have  similar  plane  expressions.  The  merging  is  done  in  two  steps.  The  first  step  is  to 
scan  each  row  of  the  point_array  looking  for  adjacent  points  in  the  same  row  which  have 
different  subregion  IDs.  For  every  pair  of  adjacent  subregions,  the  plane  expressions 
are  consulted  from  the  subregion_array  by  using  the  ID  numbers  of  the  subregions  as 
the  indices.  The  equation_difference  of  the  two  plane  equations  are  calculated  and 
compared  to  abc_difference.  If  it  less,  the  two  subregions  are  merged,  point  of  the  sec¬ 
ond  subregion  to  the  tail  point  of  the  first  subregion.  The  second  step  is  identical  except 
that  scanning  for  merging  is  done  by  column. 

D.  AN  ALGORITHM  IMPROVING  THE  CONTINUITY  OF  ADJOINING 
REGIONS 

One  common  problem  with  planar-patch  terrain  modeling  is  that  planar  patches  for 
two  adjacent  subregions  seldom  yield  a  surface  that  is  continuous  along  the  line  seg¬ 
ments  bisecting  the  points  used  to  create  the  regions.  Since  the  initial  requirements  for 
this  research  rule  out  any  other  terrain  modeling  than  planar-patch,  we  had  to  come  up 
with  some  idea  to  ensure  continuity  of  the  patches.  One  simple  solution  is  to  model 
terrain  entirely  with  small  triangles,  each  of  which  exactly  fits  three  adjacent  terrain 
points,  as  described  in  Chapter  1.  The  number  of  triangles  necessary  is  very  large. 

But  a  minor  modification  to  the  two  modeling  approaches  explored  in  this  research 
can  give  smoothness.  One  can  create  supplementary  triangular  subregions  to  fill  the  gap 
between  adjoining  subregions.  The  gap  between  adjoining  subregions  have  two  lines,  one 
for  each  side  of  the  gap.  The  lines  are  the  edges  of  the  adjacent  planes.  We  can  define 
two  triangles  on  the  gap  by  grouping  the  first  end  point  of  the  first  edge  to  the  second 
edge  and  second  point  of  the  second  edge  to  the  first  edge.  See  Figure  8  on  page  22. 
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TWO  TRIANGLES  CONNECTING 
ADJACENT  SUBREGIONS 


Figure  8.  Filling  the  Boundary  with  Two  Triangles 


This  is  better  than  the  many-trianglc  model,  but  is  still  not  optimal,  because  the  number 
of  subregions  will  be  increased  considerably  by  the  supplementary  triangular  subregions. 

Wc  have  made  a  different  attempt  to  solve  this  problem  using  smoothing  of  the  el¬ 
evation  array  and  an  expansion  to  three  dimensions  of  Rolle's  theorem.  Tne  basic  idea 
is  that  if  we  fit  planar-patches  to  an  evcrywhcre-convcx  or  everywhere-concave  region, 
the  fitted  patches  are  much  more  likely  to  be  continuous.  This  means  we  should  try  to 
ensure  that  all  the  points  within  large  regions  are  classified  as  either  planar  points,  hill 
points  or  depression  points.  See  Figure  9  on  page  23.  Smoothing  the  elevation  data 
will  help:  other  kinds  of  points  will  tend  to  be  changed  into  one  of  the  three. 

Unfortunately,  because  of  time,  we  were  not  able  to  solve  this  rather  complex 
problem  completely.  It  is  recommended  for  future  research. 


IV.  DISCUSSION 


A.  ACCURACY  ANALYSIS 

The  joint  top-down  and  bottom-up  and  strict  bottom-up  approaches  produce  fairly 
good  terrain  models  in  terms  of  accuracy  and  simplicity.  Here  accuracy  means  whether 
the  model  represented  the  real  terrain  without  losing  anything  significant,  and  simplicity 
means  the  number  of  subregions. 

The  overall  accuracy  of  a  model  can  be  verified  by  comparing  the  pictures  of  the  real 
terrain  and  model.  But  there  is  another  thing  that  describes  the  accuracy  of  a  model, 
standard  deviation.  We  can  evaluate  the  local  accuracy  in  terms  of  the  standard  devi¬ 
ation  of  the  elevation  difference  for  a  point  in  the  real  terrain  and  the  corresponding 
point  in  the  model.  In  the  top-down  approach,  each  subregion  must  have  a  fit  standard 
deviation  less  than  the  the  threshold  standard_deviation.  Let's  assume  that  the  elevation 
difference  between  a  point  on  the  real  terrain  and  a  point  on  the  model  has  a  normal 
probability  distribution.  If  we  know  the  standard  deviation  of  the  elevation  differences 
of  a  subregion: 

1.  About  68  percent  of  the  points  will  have  elevation  difference  less  than  plus  or  mi¬ 
nus  one  standard  deviation. 

2.  About  98  percent  of  the  points  will  have  elevation  difference  less  than  plus  or  mi¬ 
nus  two  standard  deviations. 

3.  Practically  all  points  (99.73  percent) 

will  have  elevation  difference  less  than  plus  or  minus  three  standard  deviations. 

So  for  example,  if  we  limit  the  model  to  have  the  worst  case  local  elevation  differ¬ 
ence,  say  10  feet,  with  9S  percent  certainty,  we  can  set  the  standard  deviation  threshold 
to  5  feet.  Consequently  every  subregion  will  have  a  standard  deviation  less  than  5  feet, 
and  98  percent  of  the  points  W'ill  have  an  elevation  difference  less  than  10  feet. 

The  threshold  standard  deviation  was  not  used  in  the  strict  bottom-up  approach  so 
the  above  argument  does  not  hold. 

B.  LIMITATIONS  OF  THE  MODEL 

In  the  joint  top-down  and  bottom-up  approach,  the  minimum  number  of  points  in 
a  subregion  is  four  and  the  shape  of  initial  subregion  must  be  a  rectangle.  This  re¬ 
striction  is  due  to  the  quadtree.  These  limitations  prohibit  the  approach  from  picking 
up  linear  (one-cell-wide)  terrain  features  and  creating  linear  subregions.  This  limitation. 
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however,  does  not  exist  in  the  strict  bottom-up  approach  owing  to  the  linked  list  of 
points  used  to  represent  a  region.  When  a  linear  region  is  a  straight  horizontal  or  ver¬ 
tical  line,  one  cannot  calculate  a  plane  equation. 

C.  EXPERIMENTAL  RESULTS 

We  now  present  the  output  of  our  terrain  modeling  program  in  graphical  form.  It  is 
applied  to  sample  elevation  data  acquired  from  the  Fort  Hunter  Liggett  terrain  data 
base.  Shown  are  the  pictures  that  are  produced  by  our  modeling  programs  using  differ¬ 
ent  thresholds. 

As  was  mentioned  in  chapter  3,  in  joint  top-down  and  bottom-up  modeling,  the 
thresholds  were  low_bound  and  upper_bound  for  gradient  magnitude,  standard_deviation, 
and  abc_difference.  In  strict  bottom-up  modeling,  the  thresholds  were  low_bound  and 
upper_b°und  of  gradient  magnitude,  curvature.  abc_difference.  and  magnitude_ratio. 

Three  sets  of  pictures  are  given  sequentially  in  the  following  section.  The  first  set 
is  produced  by  the  joint  top-down  and  bottom-up  approach  (referred  as  the  first  ap¬ 
proach  in  the  following  pictures),  the  next  by  the  strict  bottom-up  approach  (referred 
as  the  second  approach),  and  the  last  by  the  bottom-up  approach  applied  to  the 
smoothed  data  (referred  as  the  third  approach). 

In  the  array  pictures,  each  subregion  has  a  unique  code.  If  a  number  does  not  have 
a  prefix  then  the  number  is  actually  the  subregion  ID.  If  it  does  have  a  prefix  then  the 
number  has  to  be  converted.  See  the  conversion  example  below. 

23  =  =  >  subregion  ID  23. 

.23  =  =  >  subregion  ID  123. 

:23  =  =  >  subregion  ID  223. 

-23  =  =  >  subregion  ID  323. 

0  =  =  >  either  an  outermost  boundary 
or  inside  area  of  a  region. 
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Z  AXIS  (FEET) 


1.  EXPERIMENTAL  RESULTS  1  OF  THE  FIRST  APPROACH 


Figure  10.  Sample  Real  Terrain:  This  is  the  real  terrain  that  were  used  in  these 
experiments. 
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Figure  11.  Subregions  Created  by  the  Quadtree  Method  in  the  First 
Approach:  Thresholds:  lo\v_bound  =  0. 1;  upper_bound  =  0.6: 

standard_deviation=  3.  The  above  picture  shews  the  initial  subregions 
created  by  quadtree  algorithm. 
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Figure  12.  Subregions  Created  after  the  Merging  in  the  First  Approach:  Thresh¬ 
olds:  abc_difTerence=10;  standardisation  =  3;  Some  of  the  initial 
subregions  created  by  quadtree  algorithm  are  merged  by  merging 
process,  thus  yielding  less  subregions. 
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Figure  14.  Terrain  Model  Created  by  the  First  Approach  Showing  Br’jndaries 
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2.  EXPERIMENTAL  RESULTS  2  OT  FIRST  APPROACH 
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Figure  15.  Subregions  Created  by  the  Quadtree  Method  in  the  First 
Approach:  Thresholds:  Io\v_bound  =  0.1;  upper  bound  =  0.6; 

standard_dcviation=  7.  The  above  picture  shows  initial  subregions  cre¬ 
ated  by  the  quadtree  algorithm  using  a  different  threshold  (the 
standard_deviation  is  changed). 
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Figure  16.  Subregions  Created  after  the  Merging  in  the  First  Approach:  Thresh¬ 
olds:  abc_diflcrence  =  30;  standard_deviation=  7;  Some  of  the  initial 
subregions  created  by  quadtree  algorithm  are  merged  bv  the  mercing 
process.  The  thresholds  standard_de\iation  and  the  abc  dilTcrence  are 
changed. 
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Figure  18.  Terrain  Model  Created  by  the  First  Approach  Showing 
Boundaries:  The  planar  patches  are  different  from  the  planar  patches 
produced  from  the  previous  experiment. 
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3.  EXPERIMENTAL  RESl'LTS  1  OF  THE  SECOND  APPROACH 
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Figure  19.  Subregions  Created  by  the  Region  Growing  Algorithm  in  the  Second 
Approach:  Thresholds:  low_bound  =  0.1;  upper_bound  =  0.6; 

magmtude_ratio  =  0.1.  This  picture  is  produced  after  performing  the 
one-dimensional  region  growing  process  (row  region  growing). 
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Figure  20.  Subregions  Created  after  the  Merging  in  the  Second 
Approach:  Threshold:  abc_difierence  =  10;  the  two-dimensional  re¬ 
gion  growing  algorithm  is  applied  over  the  subregions  created  by  the 
one-dimensional  region  growing  algorithm  and  then  merging  is  per¬ 
formed. 
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Figure  23.  Subregions  Created  by  the  Region  Growing  Algorithm  in  the  Second 
Approach:  Thresholds:  low_bound  =  0. 1 ;  upper_bound  =  0.6; 

magnitude_ratio  =  0.2.  This  picture  is  produced  after  performing  one¬ 
dimensional  region  growing  using  a  different  threshold  (the 
magnitudc_ratio  is  changed). 
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Figure  24.  Subregions  Created  after  the  Merging  in  the  Second 
Approach:  Threshold  values  used:  abc_difTerence  =  20.  Two  dimen¬ 
sional  region  growing  is  applied  and  then  merging  is  performed  using 
different  thresholds  (abcdifferencc  and  magnitude_ration  are 
changed). 
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produced  by  the  previous  experiment. 
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Figure  28.  Subregions  Created  by  the  Third  Approach:  Thresholds: 
lo\v_bound  =  0.1;  upper_bound  =  0.6;  magnitude_ratio  =  0.1. 
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Figure  29.  Subregions  Created  after  the  Merging  in  the  Third 
Approach:  Threshold;  abc_diflercnce=  10. 
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Figure  30.  Terrain  Mode!  Created  by  the  Third  Approach:  This  picture  looks 
almost  the  same  as  the  the  picture  produced  by  the  second  approach; 
this  approach  produced  less  subregions  than  the  second  approach. 


V.  CONCLUSIONS 


In  this  research,  we  explored  three  types  of  planar-patch  terrain  modeling.  The  purpose 
of  such  models  is  a  simple  and  yet  accurate  view’  of  terrain.  Since  our  models  were  de¬ 
veloped  for  route  planning,  we  required  a  minimal  number  of  subregions,  since  this  di¬ 
rectly  affects  the  cost  of  the  route  planning. 

The  first  algorithm  (joint  top-down  and  bottom-up  algorithm)  takes  the  longest  ex¬ 
ecution  time  among  the  three.  The  second  (strict  bottom-up)  and  the  third  (bottom-up 
algorithm  applied  to  the  smoothed  elevation  data)  algorithms  take  almost  same  time, 
though  the  third  algorithm  requires  the  input  elevation  data  to  be  smoothed  first,  re¬ 
quiring  some  extra  preprocessing  time.  As  far  as  result  accuracy  (simplicity)  is  con¬ 
cerned.  the  first  algorithm  is  better  than  the  others.  The  second  and  third  do  not  use  a 
standard  deviation  threshold,  resulting  in  possible  inaccuracy,  and  tend  to  produce  more 
subregions.  Since  all  of  the  three  algorithms  use  almost  the  same  data  structure,  a  two- 
dimensional  point_array  and  a  one-dimensional  subregion_array,  they  take  almost  the 
same  space.  As  for  the  application  environments  of  these  algorithms,  the  first  algorithm 
could  be  used  in  situations  where  the  accuracy  is  critical,  while  the  second  and  third 
algoritms  could  be  preferred  in  situations  where  a  short  execution  time  is  required  and 
linear  terrain  features  should  be  identified. 

Since  the  algorithms  used  in  this  research  operate  on  large  arrays  and  involve  heavy 
computations,  a  multiple-processor  computer  architecture  would  yield  superior  process¬ 
ing  capability  since  the  same  operations  could  be  performed  for  each  part  of  the  input 
matrix  simumtancously.  As  I  conclude  this  research,  I  regret  that  I  could  not  solve  the 
continuity  issue  completely. 
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APPENDIX  A.  PASCAL  CODE  OF  JOINT  TOP-DOWN  AND 
BOTTOM-UP  ALGORITHM 


(*$s350000*) 

Author  =  Yee,  Seung  Hee. 

Date  =  3  June  1988. 

Input  =  1.  Elevation  data  acquired  from  Fort.  Hunter  Liggett 
data  base  (40  by  40  square). 

2.  file  of  magnitude  of  gradient  of  the  above  elevation 
data  created  by  using  the  program  in  Appendix  3. 

Output  =  1.  Elevation  data  of  the  model  created  by  this  program. 

2.  Elevation  data  of  the  boundary  points,  other  points 

are  given  the  value  one  so  that  when  we  display 
this  model  in  three-dimensional  pictures,  we  can 
only  see  the  boundaries  of  the  patches. 

Computer  =  Waterloo  Pascal  on  IBM  370/3033AP, 

Naval  Postgraduate  School. 

Description  = 

This  program  is  the  Pascal  implementation  of  the 
Joint  top-down  and  bottom-up  terrain  modeling. 


program  top_down( input , output ) ; 
const 

mtinft  =  3. 2808; 

low_bound  =  0.  1;  (*  mag  >  10  %  slope  then  divide  *) 

(*  mag  <=  10  %  then  disregard  *) 
upper_bound  =  0. 6;  (*  mag  <  60  %  then  do  not  subdivide  *) 

max_sd  =  3; 
abc_dif ference  =10  ; 
max_group_num=  300; 
row  =  40; 
col  =  40; 
increment  =  3; 

type 

elrange  =  1. . 10000; 
pointrec  =  record 

el  :  elrange  ; 
mag:  real; 

gp  :  0. . max_group_num; 
end; 

magrec  =  record 

mag  :  real; 
r ,c  :  1. . row; 
end; 

point_array=array  (. 1. . row, 1. . col. )  of  pointrec; 
grouppointer  =  ■■grouprec  ; 
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grouprec  =  record 

rl.rh, 

cl,cr  :  1.  .  row; 
magmax ,magmin  :  magrec  ; 
up  :  grouppointer; 
end; 

constarray  =  array  (. 1. . 3,1. . 4. )  of  real; 

nodepointer  =  ■’node; 
node  =  record 

r,c  : integer; 
next  :  nodepointer; 
end; 

coefrec  =  record 

crl,crh, 

ccl,ccr  :  1.  .  row; 

a,b,c  :  real  ; 

minmag ,maxmag  :  real; 

sd  :  real  ; 

ee,ex,ey  :  real; 

side  :  0.  .  max_group_num; 

constarr  :  constarray; 

active  :  boolean; 

list  :  nodepointer  ; 

end; 

subregion_array  =  array  ( . 1. . max_group_num. )  of  coefrec  ; 


var 

points  :  point_array  ; 
groups  :  grouprec  ; 
subregion  :  subregion_array  ; 
ok  :  constarray; 
data,datam,doutl ,dout2  :  text; 
counter , r , c , i  :  integer; 
avg,ee,ex,ey  :  real; 
top,p  :  grouppointer; 

temparr:  array  ( . 1. . max_group_num.  )  of  boolean; 

function  variance  (a,b,c,ee,ex,ey: real;  ok: constarray): real; 
begin 

variance  :=  (ee-(2*a*ex)-(2*b*ey)-(2*c*ok(.  3,4.  )) 

+( a*a*ok( .  1 , 1. ) )+( 2*a*b*ok( .2,1.)) 

+( 2*b*c*ok( .  2 , 3.  ) )+( 2*a*c*ok( .1,3.)) 
+(b*b*ok(.  2,2.  ))+(c*c*ok(. 3,3. ))) 
/(ok( .3,3. )-l); 


procedure  initialize; 
var 

r,c  : integer; 
begin 

counter  : =  0  ; 

for  c  : =  1  to  col  do 
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for  r  :  =  1  to  row  do 

points(.  r.  )(.  c.  ).  gp  :=  0  ; 
for  r  :  =  1  to  max_group_num  do 
begin 

subregion(.  r.  ).  a  :=  -999; 
subregion( . r. ) . b  :=  -999; 
subregion(.  r. ). c  :=  -999; 
subregion(. r. ). active  :=  false; 
subregion(. r. ). side  :=  0; 
subregion( .  r.  ).  sd  :=  -999; 

end; 

end;  (*  procedure  initialize  *) 


procedure  getconst(var  points: point_array;  rl, rh, cl, cr: integer; 

var  arr:  constarray;  var  ee,ex,ey  : real); 

var 

i, j: integer; 
begin 

for  i  : =  1  to  3  do 
begin 

for  j  :=  1  to  4  do 
arr(.  i, j. ): =0; 

end; 

ee  :  =  0; 
ex  : =  0; 
ey  :=  0; 

for  i  : =  rl  to  rh  do 
begin 

for  j  : =  cl  to  cr  do 
begin 

arr(. 1,1- ):=arr(. 1,1. )+( i*i); 
arr(.  1,2.  ):=arr(.  1,2. )+(j*i); 
arr(.  1,3.  ):=arr(.  1,3.  )+  i; 

arr(.  1,4.  ):  =arr(.  1,4. )+( points (. i, j. ). el*i); 

arr( .  2,1. ):=arr(. 2,1. )+(i*j); 
arr(.  2,2.  ):=arr(.2,2.  )+(j*j); 
arr( . 2 ,3.  ): =arr( . 2,3. )  +  j; 

arr(.  2,4.  ):=arr(. 2,4.  )+(points(. i, j  ). el*j); 

arr( .  3,1.  ):  =arr( .  3,1.  )+i; 
arr( .  3,2.  ):  =arr(.  3,2.  )+j; 
arr( .  3,3.  ):  =arr( .3,3.  )+l; 
arr(. 3,4. ): =arr(. 3,4. )+points(. i, j. ). el; 

ee  :=  ee  +  sqr(points(.  i, j.  ).  el); 
ex  :=  ex  +  (points( . i , j. ) . el  *  i  ); 
ey  :=  ey  +  (points(. i, j. ). el  *  j  ); 

end  (*  j  *) 

end;  (*  i  *) 

end;  (*  procedure  get  const  *) 
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procedure  gauss(ok: constarray;  var  aa,bb,cc: real); 


i,j  :  integer; 

temp  :  array  (.1.  .4.  )  of  real; 
t  :  real; 
begin 


(*  error  checking  purpose  *) 


i  •=  l; 

while  ok(.i,l.  )  B  0,0  do  i  :=  i  +  1  ;  (*  row  rotation  *) 
for  j  :=  1  to  4  do  ' 

begin 


temp(. j. )  : = 
ok(. 1, j. )  := 
ok(.  i,j. )  :  = 
end; 

t  :=  ok(.  1,4.  ); 


=  ok(.  l,j.  ); 
-  ok( .  i,  j.  ); 
=  temp(. j. ); 


for  i  :=  2  to  3  do 

for  j  :  =  4  downto  1  do 

ok(.  i,j.  )  :=  ok(. i, j.  )*ok(. 1,1. )  +  ok(. l,j. )*(-ok(. i,l. )); 
i  :  =  2; 

while  ok( . i,2. )  =  0.  0  do  i:=i+l- 

for  j  :  =  1  to  4  do 

begin 


tempi;,  j  -  ) 
ok(. 2, j. ) 
ok(. i,j.  ) 


=  ok(. .  2,  j.  ); 
=  ok(.  i,j.  ); 
=  temp(. j. ); 


for  j  : =  4  downto  1  do 

°k( .  3 ,  j .  )  :=  ok(.3,j.)*ok(.2,2.  )  +  ok(.  2,  j.  )*(-ok(.  3,2.  )); 

if  ok(. 3,3.  )  <>  0.  0  then 
begin 


ok(.  3,4.  ) 
ok(. 2,4. ) 
ok(.  1,4. ) 


aa  :=  ok(.l,4. );  (*  solution  *) 

bb  :=  ok(.2,4. );  (*  solution  *) 

cc  : -  ok(.3,4. );  (*  solution  *) 

e*dgauss  error  =  t_(  aa*ok( . 1,1. )+bb*ok(. 1,2. )+  cc*ok(. 1,3. ))  *) 
end;  (*  procedure  gauss  *) 


=  ok( . 3,4. )  /  ok( . 3,3. ); 

-  (  ok( .2,4. )  -  ok(. 3,4. )*ok(. 2,3. ))  /  ok(.2,2. ); 

-  (  ok(. 1,4. )  -  ok(. 2,4. )*ok(. 1,2. ) 

,  ^  '  °kC:3,4.)*ok(.l,3.  ))  /  okC.1,1.)  ; 


procedure  prefix(num:  integer); 
begin 

case  ((nutn  div  100)  mod  10)  of 
0  :  write( '  '); 

1  :  write( ' .  '  ); 

2  :  write( ' :  ' ); 

3  :  write( ' ); 

4  :  write( '+' ); 
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5  :  write(  '%')•, 

6  :  write( ); 

7  :  write( ); 

8  :  writ e( 1 ? ' ); 

9  :  write( '=' ); 

end; 

end;  (*  procedure  prefix  *) 

procedure  printsubl(points: point_array); 
var 

r ,c,i, j ,k  :  integer; 
begin 
page; 

for  k  : =  1  to  counter  do 
with  subregion( .  k.  )  do 

for  i  : =  crl+1  to  crh-1  do 
for  j  :=  ccl+1  to  ccr-1  do 
points(. i, j. )•  gp  :  =  0; 
for  r  :  =  1  to  row  do 
for  c  : =  1  to  col  do 
begin 

if  c<>  col  then  begin 

pref ix(points( .  r,c. ). gp); 
write( (points( . r,c. ). gp  mod  100): 2); 
end 

else  begin 

pref ix( points ( .  r,c.  ).  gp); 
writeln( (points( . r,c. ). gp  mod  100): 2) 
writeln; 
end; 

end; 

end;  (*  procedure  printsubl  *) 


procedure  printsub2(points:  point_array) ; 
var 

k,r,c,rl,r2,cl,c2  :  integer; 
current  :  nodepointer; 
begin 

for  r: =  2  to  row-1  do 
for  c: =  2  to  col-1  do 

points( . r,c. ). gp  :=  0  ; 
for  k: =  1  to  counter  do 
if  subregion(.  k. ).  active  then 

begin 

current  : =  subregion( . k. ). list  ; 
while  current  <>  nil  do  begin 
rl:  =current-.  r; 
cl:  =current-\  c; 

if  current-. next  <>  nil  then  begin 
r2: =current-. next-,  r; 
c2:  =current-.  next-,  c; 

end 

else  begin 
r2  :=  rl; 
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c2  :  =  cl; 
end; 


if 

rl 

> 

r2  then  begin 

r:  = 

r2;  r2:  = 

rl; 

rl:  —  r; 

end: 

if 

cl 

> 

c2  then  begin 

c:  = 

c2;  c2: = 

cl; 

cl:  =  c; 

end 

if 

rl! 

=r2 

then  for  c  : = 

=  cl 

to  c2  do 

points(.  rl,c.  ).  gp  :=  k; 
if  cl=c2  then  for  r  :=  rl  to  r2  do 
points(.  r,cl.  ).  gp  :=  k; 
current  :=  current-1,  next; 
end; 
end; 

for  r  : =  1  to  row  do 
for  c  :  =  1  to  col  do 
begin 

if  c<>  col  then  begin 

pref ix(points( .  r,c.  ). gp); 
write( (points( .  r,c.  ). gp  mod  100): 2); 
end 

else  begin 

pref ix(points( .  r,c.  ).  gp); 
writeln( (pointsC . r ,c. ) . gp  mod  100): 2); 
writeln; 
end; 

end; 

end; 

(*  procedure  printsub2  *) 


procedure  get_nod«s( var  subregion:  subregion_array) ; 
type 

direction  =  (  up, down, right, left  ); 

var 

i,j,k  :  integer; 

current .head ,p  :  nodepointer; 
temp  :  direction  ; 

procedure  get_next(var  p  :  nodepointer; 

var  temp  :  direction); 

var 

i,j  :  integer; 
done  :  boolean  ; 


procedure  do_dir(var  done: boolean 
begin 

new(p); 

P--  r  :=  i  ; 
p-.  c  :=  j  ; 
done  :  =  true; 
case  dir  of 

up  :  temp  : =  up; 
down  :  temp  :  =  down; 
right  :  temp  :=  right; 


dir: direction  ); 
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left  :  temp  : =  left  ; 
end;  (*  case  ,v) 

end;  (*  sub_sub_  procedure  do_dir  *) 


begin  (*  sub  procedure  get  next  *) 
dope  :=  false  ; 
i  :  =  p-’.  r; 
j  :»  P--  c; 
case  temp  of 

right  :  begin 

j  ■•= 

while  not  done  do  begin 

if  points(. i-1, j. ). gp  =  k 
then  do_dir(done,up); 
if  (not  done  and 

(points(.  i, j+1. ). gp  <>  k)) 
then  do_dir( done, down); 
if C  (subregion(. k. ). list-. r  =  i) 

and  (subregion! .  k.  ).  list-’,  c  =  j))  then 
begin 

new(p); 
p-.  r  :  =  i  ; 

P--C  :=  j  ; 
done  : =  true; 
end; 

if  not  done  then  j:=  j+1; 
end;  (*  while  *) 
end; 

left  :  begin 

j  :=  j-i; 

while  not  done  do  begin 

if  points(. i+1, j. )• gp  =  k 
then  do_dir( done, down); 
if  (not  done  and 

(points!  •  i, j -1. ).  gp  <>  k)) 
then  do_dir(done,up); 
if  not  done  then  j:=  j-1; 
end;  (*  while  *) 
end; 

up  :  begin 

i  :  =  i-1; 

while  not  done  do  begin 

if  points(.  i, j-1.  ).  gp  =  k 
then  do_dir(done,left); 
if  (not  done  and 

(points!,  i-1, j. ). gp  <>  k)) 
then  do_dir( done, right); 
if  not  done  then  i:=  i-1; 
end;  (*  while  *) 
end; 

down  :  begin 

i  :  =  i+1; 

while  not  done  do  begin 

if  points(.  i, j+1.  ).  gp  =  k 
then  do_dir( done, right); 
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if  (not  done  and 

(points( . i+1 , j. ). gp  <>  k)) 
then  do_dir( done, left); 
if  not  done  then  i: =  i+1; 
end;  (*  while  *) 
end; 

end;  (*  case  *) 

end;  (*  sub  procedure  get_next  *) 


begin  (*  procedure  get  nodes*) 
for  k  : =  1  to  counter  do 
if  subregion(. k. ). active  then 
begin 

new(p); 

p-. r  : =subregion(.  k. ). crl  ; 
p-. c  : =subregion(.  k. ).  ccl  ; 
while  (  points(.  p-*.  r-l,p-*.  c.  ).  gp  =  k)  do 

p-*.  r  :=  p-*.  r-1;  (*  move  to  the  upper  most  row  of  the  gp  *) 

while  (  points( .  p-*.  r  ,p-.  c-1.  ).  gp  =  k)  do 

p-*.  c  :=  p-’.c-l;  (*  move  to  the  left  most  col  of  the  gp  *) 

subregion( .  k.  ) .  list  :  =  p  ; 
current  :  =  p  ; 
temp  : =  right  ; 
repeat 

get_next(.  p ,  temp) ; 
current-*,  next  :=  p  ; 
current  : =  p  ; 

until  (  subregion(.  k.  ).  list-*,  r  =  p->.  r) 

and  (subregion(.  k.  ).  list-*,  c  =  p".  c)  ; 
p-*.  next  :  =  nil; 


(*  for  k  : =  1  to  counter  do 
if  subregion( . k. ). active  then 
begin 

write('gp  '  ,k:  3,  ’**=>’ ); 
current  :=  subregion( . k. ) . list  ; 
while  current  <>  nil  do  begin 

wr it e( current-*,  r:  2 , '  , '  , current-*,  c:  2 , '  =>' ); 
current  :=  current  ■*.  next; 
end; 

writeln; 

end*) 

end;(*  procedure  get_nodes  *) 


procedure  makesub(rl,rh,cl,cr  : integer;  prev: grcuppointer); 
type 

why  =  (nondivisible, magnitude, std, steep); 


var 

dif,r,c  :  integer; 
flag  :  why; 
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sd,a,b,cc  : real; 

procedure  stopdivide( f lag  :  why  ); 
var 

r,c  :  integer; 
begin 


counter  :  =  counter  +  1; 
with  subregion( .  counter. )  do 
begin 

crl  :=  rl; 
crh  :=  rh; 
ccl  : =  cl; 
ccr  : =  cr; 

maxmag  :  =  p-. magmax.  mag  ; 
minmag  :  =  p-*.  magmin.  mag  ; 
end; 

for  c  : =  cl  to  cr  do 
for  r  :=  rl  to  rh  do 

points(.  r.  )(.  c.  ). gp  :=  counter; 

(*  case  flag  of 

nondivisible  : 

write('*  non  divisible  area!'); 
magnitude  : 

write(’*  stop  divide!  by  mag.'); 
std  : 

write('*  stop,  it  fits  well.'); 
steep  : 

write('*  stop,  it  is  too  steep.'); 

end; 

write ln( ' *  area  rl,rh,cl,cr=' ,rl: 3,rh: 3, cl: 3,cr: 3, 
has  gpnum=' .counter: 3  ); 

*) 

p  :=  prev  ; 

end;  (*  sub  procedure  stop  divide  *) 


begin  (*  procedure  makesub  *) 
new(p); 

p-c  magmax.  mag  :=  -999; 

p->.  magmin.  mag  ;=  999; 

p-\  rl  :  =  rl; 

p-\  rh  :  =  rh; 

p-'.  cl  :  =  cl; 

p-'.  cr  :  =  cr; 

p-. up  : =  prev  ; 

for  c  :=  cl  to  cr  do 

for  r  :=  rl  to  rh  do 

begin 

if  points( .  r.  )( .  c.  ).  mag  >  p-c  magmax.  mag  then 
begin 

p-’.  magmax.  mag  :=  points(.  r.  )(.  c.  ).  mag; 
p-’.  magmax.  r  :  =  r; 
p->.  magmax.  c  :  =  c; 
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end; 

if  points(.  r.  )(.  c.  ).  mag  <  p-1.  magmin.  mag  then 
begin 

p**.  magmin.  mag  :=  points(.  r.  )(.  c.  ).  mag; 
p-*.  magmin.  r  :  =  r; 
p-’.  magmin.  c  :=  c; 

end; 

end; 

(*  max, min  magnitude  for  region  rl,rh,cl,cr  before  division  *) 

getconst (points , r 1 , rh, cl , cr , ok, ee , ex , ey) ; 
gauss ( ok , a , b , cc ) ; 

sd  :=  sqrt(variance(a,b,cc,ee,ex,ey,ok)); 

(*  sd  for  gp  before  division  *) 

if  ((rh  -  rl)  >  2)  and  ((cr  -  cl)  >  2)  then 
if  (p*1.  magmax.  mag  >  low_bound)  and 
(p-1.  magmin.  mag  <  upper_bound)  then 
if  sd  >  max_sd  then 
begin 

makesub(rl  ,  ((rh-rl)  div  2)  +  rl, 

cl  ,  ((cr-cl)  div  2)  +  cl,  p); 

makesub( (( rh-rl)  div  2)  +  rl  +  1  ,  rh, 
cl  ,  ((cr-cl)  div  2)  +  cl  ,  p); 

makesub( r 1  ,  ((rh-rl)  div  2)  +  rl, 

((cr-cl)  div  2)  +  cl  +  1  ,  cr  ,  p); 
makesub( (( rh-rl)  div  2)  +  rl  +  1  ,  rh, 

((cr-cl)  div  2)  +  cl  +  1  ,  cr  ,  p); 
p  :=  prev  ; 

end 

else  stopdivide(std) 
else  stopdivide( magnitude) 
else  stopdivide(nondivisible);  (*  smallest  area  *) 
end;  (*  procedure  makesub  *) 

procedure  do_combine( var  pt: point_array; 

var  subregion:  subregion_array; 

var  last, new  :  integer); 

var  ii,jj  :  integer; 

tail, current  ;  0. . max_group_num; 
begin 

subregion( . last. ). active  :=  true; 
subregion( . new. ). active  :=  false; 
current  :=  subregion(.  last.  ).  side  ; 
if  current  =  0  then 
begin 

subregion( . last. ). side  :=  new; 
current  : =  new; 

end 

else  begin 

while  current  <>  0  do 
begin 

tail  :=  current  ; 

current  :=  subregion(.  current. ). side  ; 

(*  traverse  *) 
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end;  (*  subregion  arr  * 

subregion(. tail.  ).  side  :=  new; 

current  :=  new;  (* 

end; 

while  current  <>  0  do  (*  reassigning  gp  # 
begin 

for  ii  :=  subregion( . current. ). crl  to 
subregion(. current. ).  crh  do 
for  jj  :=  subregion(.  current. ). ccl  to 
subregionC . current. ). ccr  do 
pt(. ii, jj- )• gp  :=  last; 
current  : =  subregion( . current. ). side; 
end;  (*  while  *) 
for  ii  :=  1  to  3  do 
for  j j  : =  1  to  4  do 

subregionC .  last.  ).  cons t arr ( . ii, j j. )  :  = 

subregionC- last. ). constarr(. ii,jj. )  + 
subregion(. new. ). constarr(. ii, jj. ); 
subregionC . last. ).  ee  coefs(. last. ). ee  + 

subregionC . new.  ). ee  ; 
subregionC . last. ). ex  :=  coefsC- last.  ).  ex  + 

subregionC . new. ). ex  ; 
subregionC.  last.  ). ey  :  =  coefsC- last.  ).  ey  + 

subregionC . new. )  ey  ; 
with  subregionC . last.  )  do 
begin 

gauss( const arr ,a,b,c); 

sd  :=  sqrtCvarianceCa,b,c,ee,ex,ey,constarr) ); 

end; 

end;  (*  sub  procedure  do_combine  *) 


C*  to  link  new  gp  *) 
into  existing  link  *) 

to  new  region  *) 


procedure  combine(var  pt:  point_array; 

var  subregion:  subregion_array); 

var 

tabc_difference, i , j , last , new  :  integer; 
joinabc, joinsd, joined  :  boolean; 

procedure  compareabc( var  join: boolean); 
var 

temp  :  real  ; 
begin 

temp  :=  sqrtC  sqrC subregionC . new. ). a  - 

subregionC .  last. ). a)  + 
sqr( subregionC . new. ). b  - 
subregionC . last. ). b)  + 
sqrC subregionC . new. ). c  - 
subregionC . last. ). c)  ); 
if  temp  <=  tabc_dif ference  then  join  :=  true 
else  join  : =  false; 

end;  C*  sub  procedure  compareabc  *) 

procedure  checksdC first, second: coefrec;  var  join: boolean); 
var 
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ii,jj  :  integer; 
begin 

for  ii  :=  1  to  3  do 
for  j j  : =  1  to  4  do 

first. constarr(.  )  :  = 

first. constarr( . ii,jj. )  + 
second.  constarr(. ii,jj.  ); 
first.ee  :=  first.ee  +  second.ee  ; 
first,  ex  :=  first. ex  +  second. ex  ; 
first. ey  :=  first. ey  +  second. ey  ; 
with  first  do 
begin 

gauss ( cons tarr, a, b,c); 

:  =  sqrt(variance(a,b,c,ee,ex,ey,constarr)) 

if  sd  max_sd  then  join  false  else  join  :=  true  ■ 
end;  ’ 

end;  (*  sub  procedure  check_sd_before_combine  *) 

begin  (*  procedure  combine  *) 

tabc_differer.ce  :  =  increment; 

while  t abc_d i f f er ence  <=  abc  difference  do  begin 

joined  : false; 

for  i  :=  2  to  (row-2)  do 

begin 

l.ist  :=  pt(.  i,2.  ).  gp  ; 
for  j  :=  2  to  (col-2)  do 
begin 

new  :  =  pt( .  i,  j+1.  ).  gp; 

it  ((i-row-2)  and  (j=col-2))  and  (not  joined)  then 
subregion( . new. ). active  :=  true; 
if  (newolast)  and 

( subregiori( .  new.  ) .  maxmag  <=  upper.bound)  and 
( subregion( .  last.  ).  maxmag  <=  upper_bound)  then 
begin 

compareabc( joinabc); 
if  joinabc 
then  begin 

checksd(  subregion( . last. ) , 

subregion(  .  new.  ) , joinsd  ); 
if  joinsd  then  begin 

do_combine(pt , 

subregion, last , new); 
joined  : =  true; 

end 

else  begin 

subregion( .  last.  ).  active  :=  true; 
joined  : =  false  ; 
last  : =  new  ; 

end; 

end 

else  begin 


end 


end; 


subregion( .  last. ). active  :=  true; 
joined  :=  false  ; 
last  : =  new  ; 
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else 


begin 

subregion( . last.  ). active  :=  true; 
joined  :=  false  ; 
last  :  =  new  ; 

end; 

end; 

end;  (*  row  combine  finished  *) 

for  j  : =  2  to  (row-2)  do 

begin 

last  : =  pt(.  2, j. ). gp  ; 
for  i  : =  2  to  (col-2)  do 
begin 

new  : =  pt(. i+1, j. ). gp; 
if  (new<>last)  and 

(subregion( . new. ) .  minmag  <=  upper_bound)  and 
(subregion(.  last.  ). minmag  <=  upper_bound)  then 
begin 

compareabc( joinabc); 
if  joinabc 
then  begin 

checksd(  subregion( . last. ) > 

subregion(.  new. ) , joinsd  ); 
if  joinsd  then  begin 

do_combine(pt, 

subregion, last , new); 
joined  :=  true; 

end 

else  begin 

subregion(. last.  ).  active  :=  true; 
joined  :  =  false  ; 
last  : =  new  ; 

end; 

end 

else  begin 

subregion( . last. ). active  :=  true; 
joined  :=  false  ; 
last  : =  new  ; 

end; 

end 

else  begin 

subregion( . last.  ).  active  :=  true; 
joined  : =  false  ; 
last  : =  new  ; 

end; 

end; 

end;  (*  column  combine  finished  *) 
tabc_dif ference  :=  tabc_dif ference  +  increment; 
end;  (*  while  *) 

end;  (*  procedure  combine  *) 

procedure  cal_avg_dev( var  subregion: subregion_array; 

var  points: point_array;  var  avg:real); 

var 

sum:  real; 
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1, j ,gp: integer; 
begin 

sum  :  =  0.0; 

for  i: =  2  to  row-1  do 
for  j:=  2  to  col-1  do 
begin 

gp  : =  points(. i, j. ). gp; 

sum  :=  sum  +  abs((  subregion( . gp. ).  a*i  + 
subregion( . gp. ) . b*j  + 
subregion(. gp. ). c) 

-  points(. i, j. ). el); 

end; 

avg  :=  sum  /  ( ( row-2)*( col-2) ) ; 
writeln( ' average  deviation  =  ' ,avg  : 15); 
end;  (*  sub  procedure  cal_average_deviation  *) 


procedure  writedata(points: point_array); 
var 

r,c  :  integer; 
tempa , tempb , tempc  :  real; 
begin 

for  r: =  2  to  row-1  do 
for  c: =  2  to  col-1  do 
begin 

tempa  :=  subregion( .  points( .  r , c. ) . gp. ) . a; 
tempb  :=  subregion( .  points( . r , c. ) . gp. ) . b; 
tempc  :=  subregion( . points( . r , c. ) . gp. ) . c; 
points(. r,c. ). el  :=  trunc  (tempa*r  +  tempb*c  +  tempc); 
end; 

for  r,  :  =  1  to  col  do 
for  r  : =  row  downto  1  do 

writeln(doutl ,points( . r,c. ). el); 

end; 

(*  procedure  writedata  *) 


procedure  writetempdata(points:  point_array ) ; 
var 

k , r , c , rl , r2 , cl ,c2  :  integer; 
tempa , tempb , tempc  :  real; 
current  :  nodepointer; 
begin 

for  r: =  2  to  row-1  do 
for  c: =  2  to  col-1  do 

points(. r,c. ). el  :=  1  ; 
for  k: =  1  to  counter  do 
if  subregion(. k. ). active  then 

begin 

current  :=  subregion(. k. ). list  ; 
while  current  <>  nil  do  begin 
rl:  =current-.  r; 
cl:  =current-’.  c; 

if  current-1,  next  <>  nil  then  begin 
r2:  =current-’.  next-1,  r; 
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c2:  =current-'.  next-*,  c; 

end 

else  begin 
r2  :=  rl; 
c2  :  =  cl; 
end; 

if  rl  >  r2  then  begin  r: =  r2;  r2: =  rl;  rl: =  r;  end; 


if  cl  >  c2  then  begin  c:  =  c2;  c2:=  cl;  cl:  -  c;  end; 

teropa  :=  subregion( . k. ). a; 
tempb  :=  subregionC . k. ) . b; 
tempc  :=  subregion( . k. ) . c; 
if  rl=r2  then  for  c  :  =  cl  to  c2  do 

points(. rl,c. ). el  :=  trunc  (tempa*rl  +  tempb*c  +  tempc); 
if  cl=c2  then  for  r  :=  rl  to  r2  do 

points(. r,cl. ). el  :=  trunc  (tempa*r  +  tempb*cl  +  tempc); 
current  :=  current ■*.  next; 
end; 
end; 

for  c  : =  1  to  col  do 
for  r  : =  row  downto  1  do 

writeln(dout2 ,points( . r,c. ). el); 

end; 

(*  procedure  writetempdata  *) 


procedure  check; 
var 

num,r,c,i  :  integer; 
begin 

for  c:  =  1  to  counter  do  temparr(.c. )  :=  false;  (*  check  *) 
num  : =  0; 

for  i  :  =  1  to  counter  do 

if  subregion( . i.  ).  active  then 
begin 

r:  =  i  ; 

temparr(.r.  )  :=  true  ; 
num  : =  num  + 1  ; 
while  r  <>  0  do 
begin 

with  subregionC . r. )  do 
temparr( . r. )  :=  true; 
r: =  subregionC • r. ). side  ; 

end; 

end; 

for  c:  =  1  to  counter  do  temparrC.c. )  : =  false;  C*  check  *) 
num  :=  0; 

for  i  : =  1  to  counter  do 

if  subregionC. i.  )•  active  then 
begin 

r:=  i  ; 

temparrC-r.  )  :=  true  ; 
num  : =  num  +1  ; 
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while  r  <>  0  do 
begin 

with  subregion( .  r. )  do 
temparr( . r.  )  :  =  true; 
r: =  subregion(. r. )• side  ; 

end; 

end; 


writeln('num  of  gps  =  ’ , 
1  upper_bound= 
writeln( ’ unexamined  gps 
c:=  0; 


counter  :3,’  *  lowbound=' , low_bound: 3, 
, upper  bound:  3, 1  *  max_sd='Tniax  sd:  3 

are  =>T); 


for  i  : =  1  to  counter  do  if  temper r( . i. )  =  true  then  c: =  c+1 
else  write('*' ,i:3, 
writeln; 

of  8ps  examined  =>’,0:4,’  *  abc  difference*'  , 
abc_dif ference); 

writeln('#  of  combined  gps  =>' ,num  :4); 
end;  (*  procedure  check  *) 
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Main  Program 


(* 


===== - - - ==========================*) 

begin 

page; 

initialize  ; 

new(top); 

top-1,  up  :  =  nil  ; 

reset(data, 'data  ell'); 
reset(datam, 1  data  magi'); 

rewrite(doutl, 'data  outl');  (*  data  after  combine  process  *) 
rewrite(dout2, ' tdata  outl');  (*  data  of  boundary  *) 

for  c  : =  1  to  col  do 
for  r  : =  row  downto  1  do 

readln(data,points(. r, )(. c. ). el); 
for  c  : =  1  to  col  do 
for  r  : =  row  downto  1  do 

readln(datam,points( .  r.  )(.  c.  ).  mag); 

makesub(2,(row-l) ,2, (col-1), top); 

printsubl(points);  (*  print  the  subregions  created  by  makesub  *) 

for  i  : =  1  to  councer  do 
with  subregion(. i. )  do 
begin 

getconst( points , cr 1 ,crh, ccl ,ccr ,ok,ee ,ex,ey); 
gauss(ok,a,b,c); 

sd  :=  sqrt(variance(a,b,c,ee,ex,ey,ok)); 

(*  sd  for  region  i  after  division  *) 

constarr  :  =  ok  ; 
end; 

combine( points .subregion); 
get_nodes( subregion) ; 

printsub2(points);  (*  print  the  subregions  created  by  combine  *) 
writedata(points)  ; 

(*  write  modified  elevation  data  into  file  'data  out#'*) 
wr it etempdata( points ) ; 

(*  write  elevation  data  into  file  'tempdata  out#'  *) 
cal_avg_dev( sub region, points , avg); 

(*  check;  *) 

end. 
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APPENDIX  B.  PASCAL  CODE  OF  STRICT  BOTTOM-UP  ALGORITHM 


(*$s450000*) 

( ************************ 


Author 

Date 

Input 


Output 


Computer  = 


Yee,  Seung  Hee. 

3  June  1988. 

1.  Elevation  data  acquired  from  Fort.  Hunter  Liggett 
data  base  (40  by  40  square). 

2.  file  of  magnitude  of  gradient  of  the  above  elevation 
data  created  by  using  the  program  in  Appendix  3. 

3.  file  of  curvature  of  points  created  by  using  the 
program  in  Appendix  3. 

1.  Elevation  data  of  the  model  created  by  this  program. 

2.  Elevation  data  of  the  boundary  points,  other  points 
are  given  the  value  one  so  that  when  we  display 
this  model  in  three-dimensional  pictures,  we  can 
only  see  the  boundaries  of  the  patches. 

Waterloo  Pascal  on  IBM  370/3033AP, 

Naval  Postgraduate  School. 


Description  = 


This  program  is  the  Pascal  implementation  of  the 
Joint  top-down  and  bottom-up  terrain  modeling. 


"fcirki c  ) 


PROGRAM  Bottom_up( input .output); 
const 

low_bound  =  0. 1;  (*  magnitude  <  10  %  then  level  slope  *) 

(*  magnitude  >  50  %  then  unsafe  slope  *) 

up_bound  =  0.6;  (*  else  safe  slope  *) 

magnitude_ratio  =  0.2;  (*  ratio  of  magnitude  : 

dif  of  two  mags  /  old  mag  *) 

abc_comparison_const  =20; 
max_group_num  =  400; 
row  =  40; 
col  =  40; 
type 

elrange  =  1. . 10000; 

mag_type  =  ( level, safe, unsafe); 

point  =  record 

r,c: integer; 
end; 

nodepointer  =  ■•node; 
node  =  record 

r,c  : integer; 
next  :  nodepointer; 
end; 

point_rec  =  record 

el  :  elrange  ; 
mag  :  real; 
magtype:  mag  type: 
cur:  char; 

gp  :  0. .  max_group_num; 
ptnext  :  point; 
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end; 

point_array  =  array  (.  1.  .  row, 1. . col. )  of  point_rec; 
polynomial_array  =  array  (. 1. . 3, 1. . 4. )  of  real; 
coef_record  =  record 

start  :  point; 

upper_left  :  point; 

a,b,c  :  real  ; 

sd  :  real  ; 

ee,ex,ey  :  real; 

side  :  0. . max_group_num; 

poly_array  :  polynomial_array; 

active  :  boolean; 

list  :  nodepointer  ; 

numpts  :  integer; 

avgmag: real; 

end; 

coe^— array  —  array  (. 1. . max_group_num. )  of  coef  record  ; 
var  - 

points  :  point_array  ; 
subregion  :  coef_array  ; 
ok  :  polynomial_array; 

data,datain,datac,doutl,dout2,dout3Jdout4  :  text; 
counter, r ,c, i,num  :  integer; 
tempmag  :  real; 

function  variance  (a,b,c,ee,ex,ey: real;  ok: polynomial  array):real- 
begin  ~ 

variance  :=  ( ee-( 2*a*ex) -( 2*b*ey) -( 2*c*ok( . 3 ,4. ) ) 

+( a*a*ok( .  1 , 1.  ) )+( 2*a*b*ok(  .2,1.)) 
+(2*b*c*ok(. 2,3.  ))+(2*a*c*ok(. 1,3. )) 
+(b*b*ok(.  2,2.  ))+(c*c*ok(. 3,3. ))) 

/ ( ok( . 3,3. ) *  1) ; 

end; 

procedure  initialize; 
var 

r ,c  : integer; 
begin 

counter  : =  0  ; 
for  c  :=  1  to  col  do 
for  r  : =  1  to  row  do 
begin 

points(.  r , c. ) . gp  :=  0  ; 
points( .  r ,c.  ).  ptnext. r  :  =  0  ; 
points( .  r,c.  ).  ptnext.  c  :=  0  ; 
end; 

for  r  : =  1  to  max_group_num  do 
begin 

subregion( .  r.  ).  a  :=  -999; 
subregion( .  r.  ).  b  :=  -999; 
subregion( .  r.  ).  c  :=  -999; 
subregion(. r.  ).  active  :=  true; 
subregion(.  r. ).  side  :=  0; 
subregion(.  r.  ).  sd  :=  -999; 


end; 


(*  procedure  initialize  *) 
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procedure  getconst(  r,c: integer; 

var  arr: polynomial_array;  var  ee,ex,ey  : real) 

begin 

arr(.  1,1. ):  =arr(. 1,1. )+(r*r); 
arr( .  1,2.  ):=arr(.  1,2.  )+(c*r); 
arr(. 1,3. ):=arr(. 1,3. )+  r; 

arr( .  1,4.  ):  =arr( .1,4. )+( points ( . r ,c. ). el*r); 

arr(.  2,1.  ):=arr(.  2,1.  )+(r*c); 
arr(.  2,2.  ):  =arr(.  2,2. )+(c*c); 
arr( .  2 ,3. ): =arr( . 2 , 3. )  +  c; 

arr(.  2,4.  ):=arr(.  2,4.  )+(points(. r,c. ). el*c); 

arr( .  3,1.  ):  =arr(  .3,1.  )+r; 
arr( .3,2.):  =arr( .3,2.  )+c; 
arr(.  3,3.  ):  =arr(.  3,3.  )  +  '; 
arr( .3,4.):  =arr( .3,4.  )+points( .  r ,c.  ).  el; 

ee  : =  ee  +  sqr(points(. r,c. ). el); 
ex  :  =  ex  +  ( points ( . r,c. ). el  *  r  ); 
ey  :=  ey  +  (points(. r,c. ). el  *  c  ); 
end;  (*  procedure  get  const  *) 

procedure  gauss(ok:  polynomial_array;  var  aa,bb,cc:  real); 


var 

i,j  :  integer; 

temp  :  array  (.1.  .4.  )  of  real; 

t  :  real;  (*  error  checking  purpose  *) 

begin 

i  :=  1; 

while  ok( . i , 1. )  =  0. 0  do  i  :=  i  +  1  ;  (*  row  rotation  *) 

for  j  : =  1  to  4  do 

begin 

temp(. j. )  : =  ok( . 1 , j. ); 
ok(.  1, j.  )  :  =  ok(.  i, j.  ); 
ok( . i, j. )  : =  temp( . j. ); 
end; 

t  :  =  ok(. 1,4. ); 

for  i  :  =  2  to  3  do 

for  j  : =  4  downto  1  do 

ok(. i,j. )  :=  ok( . i , j . )*ok(. 1,1. )  +  ok(. l,j. )*( -ok( . i , 1. ) ); 
i  :  =  2; 

while  ok(.i,2.  )  =0.0  do  i:=i  +  l; 

for  j  : =  1  to  4  do 

begin 

temp( . j. )  : =  ok(. 2, j. ); 
ok(. 2, j. )  : =  ok( . i, j. ); 
ok( . i, j. )  :=  tempC . j. ); 
end; 

for  j  : =  4  downto  1  do 
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ok(.  3 ,  j .  )  :  -  .  k(. 3,j.  )*ok(.  2,2.  )  +  ok(.  2,j.  )*(-ok(.  3,2.  )); 


if  ok(. 3,3. )  <>  0.  0  then 
begin 


ok(. 3,4. ) 
ok(. 2,4. ) 
ok( . 1,4. ) 


aa  :  =  ok(. 1,4. ) 
bb  :=  ok(. 2,4. ) 
cc  : =  ok(.  3,4.  ) 


=  ok( .3,4.  )  /  ok(.  3,3.  ); 

-  (  ok( .2,4. )  -  ok(. 3,4. )*ok(. 2,3. ))  /  ok(.2,2.  ) 
=  (  ok(. 1,4.  )  -  ok(. 2,4.  )*ok(.  1,2.  ) 

-  ok(.  3,4.  )*ok(.  1,3.  ))  /  okC.1,1.)  ; 


(*  solution  *) 

(*  solution  *) 

.  -  - ,  -  , ,  (*  solution  *) 

(*writeln( 1  gauss  error  =' , 

t-(  aa*ok(. 1,1.  )+bb*ok(.  1,2. )+  cc*ok(. 1,3. )))*) 

end; 

end;  (*  procedure  gauss  *) 


procedure  prefix(num: integer); 
begin 

case  ((num  div  100)  mod  10)  of 
0  :  write( '  ' ); 

1  :  write( ' . 1 ); 

2  :  write( ' : ' ); 

3  i  write( 1  - ' ); 

4  :  write( 1 +'); 

5  :  write( ’ %' ); 

6  :  write('*'); 

7  :  write( ' &' ); 

8  :  write( ' ? ' ); 

9  :  write( ' =' ); 

end; 

end;  (*  procedure  prefix  *) 

procedure  printsubcpoints: point_array); 
var 

r,c  :  integer; 
begin 
page; 

for  r  : =  1  to  row  do 
for  c  : =  1  to  col  do 
begin 

if  c<>  col  then  begin 

pref ix(points( . r,c. ). gp); 
write((points(. r,c. ). gp  mod  100): 2); 
end 

else  begin 

pref ix(points( .  r,c. ). gp); 
wr it e ln( (point s( . r ,c. ) . gp  mod  100): 2); 
writeln; 
end; 

end; 

end;  (*  procedure  printsub  *) 


procedure  printsub2(points: point_array); 
var 

k,r ,c,rl ,r2,cl,c2  :  integer; 
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current  :  nodepointer; 
begin 
page; 

for  r: =  2  to  row-1  do 
for  c: ~  2  to  col-1  do 

points( .  r , c.  ).  gp  :=  0  ; 
for  k: =  1  to  counter  do 
if  subregion( . k.  ). active  then 
begin 

current  :=  subregion( . k.  ).  list  ; 
while  current  <>  nil  do 
begin 

rl:  =currenf.  r; 
cl:  ^current-*.  c; 

if  current-*,  next  <>  nil  then  begin 
r2:  =currenf.  next-1,  r; 
c2:  =currenf.  next-1,  c; 

end 

else  begin 
r2  :=  rl; 
c2  :=  cl; 
end; 

if  rl  >  r2  then  begin  r:  =  r2;  r2: =  rl;  rl: =  r;  end; 
if  cl  >  c2  then  begin  c: =  c2;  c2: =  cl;  cl: =  c;  end; 

if  rl=r2  then  for  c  :=  cl  to  c2  do 
points(. rl,c.  ).  gp  :  =  k; 
if  cl=c2  then  for  r  : =  rl  to  r2  do 
points( .  r , cl.  ) .  gp  :=  k; 
current  :=  current-*,  next; 
end; 
end; 

for  r  : =  1  to  row  do 
for  c  : =  1  to  col  do 
begin 

if  c<>  col  then  begin 

prefix(points(. r,c. ). gp); 
write( (points( . r,c. ). gp  mod  100): 2); 
end 

else  begin 

prefix(points( . r,c. ). gp); 
writeln( (points( . r,c. ). gp  mod  100): 2); 
wr iteln; 
end; 

end; 

end; 

(*  procedure  printsub2  *) 

procedure  printmagcur; 
var 

r,c  :  integer; 
begin 
(*  page*) 

for  r  : =  1  to  row  do 
for  c  :  =  1  to  col  do 
begin 


70 


end; 


if  c<>  col  then  begin 

case  points(.  r,c.  ).  magtype  of 

level  : write( ' 1 ' ,points(.  r ,c.  ).  cur , '  1 
safe  : write( ' s j ,points( .  r , c.  ).  cur , '  ' 
unsafe  : write( ' u*  , points ( . r , c.  ).  cur, ' 
end;  (*  case  *) 
end 

else  begin 

case  points(.  r,c.  ).  magtype  of 

level  : writeln( ' 1' , point s( .  r ,c.  ).  cur, ' 
safe  : write ln( ' s ' ,points( . r ,c. ). cur, 1 
unsafe: writeln( 'u' ,points(.  r,c.  ).  cur, ' 
end; 

writeln; 

end; 

end; 

(*  procedure  printmagcur  *) 


); 

); 

j 


); 


); 

); 

); 


procedure  makeone_rowl(var  pointl ,point2  :  point); 
begin 

points(. pointl. r, pointl. c. ). ptnext. r  point2. r; 
points( .  pointl.  r, pointl. c.  ).  ptnext. c  :  =  point2. c; 
points( . point2.  r,point2.  c.  ).  gp  :  = 

point s( .  pointl.  r , pointl.  c.  ).  gp; 
with  subregion( . counter. )  do 
numpts  : =  numpts  +1; 


end;  (*  sub  procedure  make  one_l  *) 


procedure  makeone_row2( var  pointl ,point2  :  point); 
begin 

points( .  pointl.  r, pointl.  c.  ).  ptnext.  r  :=  pointZ.  r; 
points( .  pointl.  r, pointl.  c.  ).  ptnext.  c  :  =  point2. c; 
points( . point2.  r, pointZ.  c.  ).  gp  :  = 

points(.  pointl.  r  .pointl.  c.  ).  gp; 
with  subregion( . counter.  )  do 
begin 

numpts  :=  numpts  +1; 

avgmag  :=  avgmag*( numpts- 1) /numpts  + 

points(.  point2.  r,point2.  c. ). mag/numpts; 

end; 

end;  (*  sub  procedure  make  one_2  *) 


procedure  rowlink; 
var 

i, j  :  integer; 

pointl ,point2  :  point; 
begin  (*  sub  procedure  rowlink  *) 

counter  :=  0; 
for  i  : =  2  to  (row-1)  do 
begin 
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counter  :=  counter  +  1  ; 
points (.  i, 2.  ).  gp  : =  counter; 
subregionC- counter. ). start. r  :=  i; 
subregionC. counter. ).  start. c  :=  2; 
subregion( . counter. ) . upper_lef t  : = 

subregionC . counter. ). start  ; 
subregion(.  counter. ).  avgmag  :=  points (. i, 2. ). mag; 
subregionC . counter.  ).  numpts  :=  1  ; 
for  j  :=  2  to  C col-2)  do 
begin 

pointl. r  :  =  i; 
point  1. c  :  =  j; 
point2.  r  :  =  i; 
point2.  c  : =  j+1; 

if  CCpointsC- i, j+1. )•  magtype  <>  safe)  and 

(pointsC. i, j. )■ magtype  =  pointsC- i, j+1. ). magtype) )  then 
makeone_row 1 C  po int 1 , point  2 ) 

(*  C level  level)  or  Cunsafe  unsafe)  *) 
else  if  ((pointsC.  i,j. ).  magtype  =  safe)  and 

(pointsC . i , j+1. ). magtype  =  safe))  and 
CpointsC- i, j-  ).  cur  =  pointsC-  i, j+1.  ).  cur)  then 
if  (abs(points(. i, j+1.  )■  mag- 

subregionC-  counter.  ).  avgmag)/ 
subregion(..  counter.  ).  avgmag)  < 
magnitude_ratio  then 
makeone_row2( pointl ,point2) 
else  begin 

counter  : =  counter  +1  ; 
points( .  i, j  +  1.  ).  gp  :=  counter; 
subregionC .  counter. ). start  :=  point2; 
subregionC .  counter.  ). upper_le ft  :=  point2  ; 
subregionC . counter. ) .  avgmag  :  = 

pointsC- i, j+1.  )•  mag; 
subregionC . counter. ). numpts  :=  1  ; 
pointl  :=  point2; 

end 

else  begin 

counter  :  =  counter  +1  ; 
pointsC- i, j+1.  )•  gp  :  =  counter; 
subregionC .  counter.  ).  start  :=  point2; 
subregionC- counter. ). upper_left  :=  point2  ; 
subregionC . counter.  ) .  avgmag  :  = 

pointsC- i, j+1.  )•  mag; 
subregionC . counter. ). numpts  :=  1  ; 
pointl  :=  point2; 

end; 

end; 

end; 

end;  (*  sub  procedure  row  link  *) 


procedure  makeone_colC last , new  :  point); 
var 

ii,jj  :  integer; 
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tail, current  :  point; 


begin 

subregion(.  points(.  new.  r, new.  c.  ).  gp.  ).  active  :=  false  ; 
with  subregion( . points( . last. r, last.  c. ). gp. )  do 
begin 

numpts: =numpts  + 

subregion(.  points(.  new.  r.new.  c.  ). gp.  ).  numpts; 
avgmag:  =avgmag*( numpts - 

subregion(.  points (.  new.  r.new.  c.  ). gp.  ).  numpts)/ 
numpts  + 

subregion(. points(. new. r.new. c. ). gp. ). avgmag* 
subregion(.  points(. new. r.new. c. ). gp.  ).  numpts /numpts; 
end;  (*  reassigning  number  of  points  and  average  magnitude  *) 

if  (subregion(.  points(.  new.  r.new.  c. ). gp.  ). upper_left. r  + 
subregion(.  points(.  new.  r.new. c. ). gp.  ). upper_left. c)  < 

(subregion(  .  point s( .  last,  r , last. c. ).  gp. ). upper_left. r  + 
subregior( .  point s( .  last,  r.last. c. ). gp. ). upper_left. c)  then 
subregion(.  points(.  last. r, last. c.  ). gp.  ). upper_left  :  = 
subregion( .  points( . new. r.new. c. ). gp.  ). upper_left  ; 

current  :=  points(  .  last,  r , last. c. ). ptnext  ; 

if  current. r  =  0  then  (*  or  current. c  =  0  *) 

begin 

points ( .  last,  r, last.  c.  ).  ptnext  :  = 

subregion( . points( . new. r.new. c. ) . gp. ) . start  ; 
current  :=  subregion( .  points( . new. r .new. c.  ). gp.  ). start  ; 

end 

else  begin 

while  current. r  <>  0  do 
begin 

tail  : =  current; 

current  :=  points( . current. r , current,  c.  ).  ptnext; 
end;  (*  traverse  ptnext  of  lastpt  link 

to  the  end  of  link  *) 

points(. tail. r, tail. c.  ).  ptnext:  =  (*  now  link  the  tail 

points  to  the  head  of*) 
subregion( .  point s( . new. r.new. c. ).gp. ). start  ; 

(*  newpt  link  *) 

current  :=  (*  into  existing  link  *) 

subregion(.  points(. new. r.new. c. ). gp. ). start  ; 

end; 

while  current. r  <>  0  do  (*  reassigning  gp  #  to  new  pts  *) 
begin 

points(. current. r, current. c. ). gp  :=  points( . last. r, last. c. ). gp; 
current  :=  points(.  current. r, current. c. ). ptnext; 
end;  (*  while  *) 

end;  (*  procedure  makeone_col  *) 


procedure  collink; 
var 

i.j  :  integer; 
pointl ,point2  :  point; 
begin  (*  sub  procedure  col  link  *) 
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for  j  : =  2  to  (col-1)  do 
begin 

for  i  :  =  2  to  (row-2)  do 
begin 

pointl.  r  :  =  i; 
point  1.  c  :  =  j ; 
point2. r  : =  i+1; 
point2. c  : =  j ; 

if  (points(. i+1, j. ). gp  <>  points( . i , j . ) . gp)  then 
if  ( ( points( . i+1 , j . ) . magtype  <>  safe)  and 

(points(. i, j. ). magtype  =  points! . i+1 , j .)• magtype) )  then 
makeone_col(pointl ,point2) 

(*  ( level , level)  or  ( unsafe , unsafe)  *) 
else  if  ( (points! • i , j- )■ magtype  =  safe)  and 

(points!  •  i+1 ,  j- )• magtype  =  safe))  then 
if  (points(. i, j. ). cur  =  points( . i+1 , j . ) . cur )  then 
if  abs( subregion( .  points ( . i+1 , j .  ) . gp.  ) . avgmag - 
subregion! . points( . i , j. ) . gp. ) . avgmag) / 
subr egion( .  points(  .  i , j .  ) .  gp. ) . avgmag  < 
magnitude_ratio 

then 

makeone_col( pointl ,point2); 
pointl  :  =  point2; 
end; 
end; 

end;  (*  sub  procedure  collink  ,v) 
procedure  do_combine( last , new  :  point); 


ii.jj  :  integer; 
tail, current  :  point; 


begin 

subregion! .  points( .  new. r , new. c.  ). gp.  ). active  :=  false  ; 
with  subregion( .  points( .  last,  r , last. c. ). gp. )  do 
begin 

numpts:  =numpts  + 

subregion!,  points!,  new. r ,new. c.  ). gp.  ). numpts; 
end;  (*  reassigning  number  of  points  *) 

if  ( subregion! .  points! .  new.  r , new.  c.  ).  gp.  ).  upper_left.  r  + 
subregion! . points! .  new.  r ,new.  c.  ).  gp.  ) .  upper_left.  c)  < 
(subregion! .  points!  •  last. r , last. c. ). gp. ). upper_left.  r  + 
subregion! .  points! •  last.  r,last.c. ).gp. ). upper_left. c)  then 
subregion!  •  points! .  last,  r , last.  c.  ) . gp. ) .  upper_left  :  = 
subregion!  .points! .  new.  r , new.  c.  ) .  gp.  ) .  upper_lef t  ; 


current  :=  points! . last. r , last. c. ). ptnext  ; 

if  current,  r  =  0  then  (*  or  current. c  =  0  *) 

begin 

points!. last. r , last.  c.  ).  ptnext  :  = 

subregion! . points! . new. r.new.c. ).gp. ). start  ; 
current  :=  subregion! . points! .  new. r , new. c. ). gp.  ). start  ; 


end 


else  begin 

while  current. r  <>  0  do 
begin 

tail  : =  current; 

current  :=  points( . current . r , current. c. ). ptnext; 
end;  (*  traverse  ptnext  to  link  the  head  *) 

points(  .  tail,  r , tail.  c. ). ptnext: =  (*  of  the  *) 

subregion(.  points(. new. r.new. c. ). gp.  ). start  ; 

(*  new  point  ») 

current  :=  (*  into  existing  link  *) 

subregion(. points(.  new. r,new. c. '. gp. ). start  ; 

end; 

while  current. r  <>  0  do  (*  reassigning  gp  #  to  new  pts  *) 
begin 

points( .  current,  r .current. c. ). gp  :=  points( .  last,  r , last. c. ). gp; 
current  :=  points( . current. r , current. c. ). ptnext; 
end;  (*  while  *) 


for  ii  : =  1  to  3  do 
for  j j  : =  1  to  4  do 

subregion( . points( .  last,  r , last.  c. ).  gp 
subregion( .  points (  .  last,  r , last,  c 
subregion( . points( .  new.  r .new.  c.  ) 
subregion( . point s( .  last,  r , last.  c.  ).  gp 
subregion( .  point s( .  last,  r , last,  c 
subregion( . points( .  new.  r .new.  c.  ) 
subregion( . points ( . last,  r , last.  c. ) .  gp 
subregion( .  points ( .  last. r , last,  c 
subregion( .  points( .  new.  r .new.  c.  ) 
subregion(. points(. last,  r , last.  c.  ).  gp 
subregion( . point s( . last,  r , last,  c 
subregion( . point s( .  new.  r .new.  c.  ) 
with  subregion( . points ( .  last,  r , last,  c 
begin 


). poly_array(.  ii, j j.  )  :  = 

). gp.  ). poly_array(.  ii, j j.  )  + 
gp. ). poly_array(.  ii, j j-  ); 

)  .  ee  :  — 

).  gp.  ).  ee  + 
gp.  ). ee  ; 

) .  ex  :  = 

).  gp.  ).  ex  + 
gp.  ).  ex  ; 

) .  ey  :  = 

).  gp.  ).  ey  + 
gp-  )• ey  ; 

).  gp.  )  do 


sd  :=  sqrt(variance(a,b,c,ee,ex,ey,poly_array)); 


end; 

end;  (*  procedure  do_combine  *) 


procedure  combine(var  pt:  point_array; 

var  subregion: coef_array) ; 

var 

i.j  :  integer; 

last, new  :  point; 

joinabc ,  joinsd , joined  :  boolean; 

procedure  compareabc( var  join:  boolean); 
var 

temp  :  real  ; 
begin 

temp  :=  sqrt(  sqr(subregion( .  poir.ts(.  new.  r  .new.  c.  ).  gp.  ).  a  - 
subregion( .points(. last. r.last. c. ). gp. ). a)  + 
sqr( subregion( . points ( . new. r .new. c. ) . gp. ) .  b  - 
subregion( . points( . last. r, last. c. ). gp. ). b)  + 
sqr(  subregior.(  .  points( .  new.  r  .new.  c.  ).  gp.  ) .  c  - 


subregion( . points(. last. r.last.  c. ).  gp.  ).  c)  ); 
if  temp  <=  abc_compar ison_const  then  join  :=  true 
else  join  :=  false; 

end;  (*  sub  procedure  compareabc  *) 


begin  (*  procedure  combine  *) 

joined  :=  false; 
for  i  : =  2  to  (row-2)  do 
begin 

last. r  :  =  i; 
last. c  :  =  2; 

for  j  :  =  2  to  (col-2)  do 
begin 

new.  r  :  =  i; 
new.  c  : =  j+1; 

if  ( ( i=row-2)  and  (j=col-2))  and  (not  joined)  then 

subregion( . points( . new. r ,new. c. ) . gp. ) . active  : = 

true; 

if  ( points( . new. r ,new. c. ) . gp  <>  points( . last. r , last. c. ). gp) 
then  begin 

compareabc( joinabc); 
if  joinabc 
then  begin 

do_combine( last, new); 
joined  : =  true; 
end 

else  begin 

subregion( .  point s( .  last. r , last. c. ). gp. ). active  :  = 

true; 

joined  : =  false  ; 
last  :=  new  ; 
end;  (*  joinabc  *) 

end 

else  begin 

subregion( . points(. last. r.last. c. ).  gp. ).  active  :  = 

true; 


joined  :  =  false  ; 
last  : =  new  ; 

end;  (*  last  <>  new  *) 
end; (*  end  j  *) 

end;  (*  end  i  *)(*  row  combine  finished  *) 

for  j  : =  2  to  (row-2)  do 

begin 

last. r  : =  2; 
last,  c  : =  j; 

for  i  :=  2  to  (col-2)  do 
begin 

new. r  : =  i+1; 
new.  c  : =  j  ; 

if  (points(. new. r ,new. c. ). gp  <>  points( . last. r , last. c. ). gp) 
then  begin 

compareabc( joinabc); 
if  joinabc 
then  begin 
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do_combine( last , new) ; 
joined  : =  true; 

end 

else  begin 

subregion(  .  points( . last,  r .last.  c.  ).  gp.  ).  active 

true; 

joined  :  —  false  ; 
last  : =  new  ; 
end;  (*  join  abc  *) 

end 

else  begin 

subregion(. points(.  last,  r, last.  c.  ).  gp.  ).  active  :  = 

true; 

joined  :=  false  ; 
last  : =  new  ; 

end;  (*  last  <>  new  *) 

end; 

end;  (*  column  combine  finished  *) 
end;  (*  procedure  combine  *) 


procedure  get_boundaries( var  subregion:  coef  array)- 
type 

direction  =  (  up .down , right , left  ); 

var 

i. j ,k,r ,c  :  integer; 

current .head, p  :  nodepointer; 
temp  :  direction  ; 

procedure  get_next(var  p  :  nodepointer; 

var  temp  :  direction); 

var 

i.j  :  integer; 
done  :  boolean  ; 


procedure  do_dir(var  done:  boolean  ;  dir: direction  ); 
begin 

new(p); 
p-.  r  :  =  i  ; 

P".  c  :  =  j  ; 


done  :  = 

true; 

case  dir 

of 

up 

:  temp 

=  up; 

down 

:  temp 

=  down; 

right 

:  temp 

=  right; 

left 

:  temp 

=  left  ; 

end;  (*  case  *) 

end;  (*  sub_sub_  procedure  do_dir  *) 


begin  (*  sub  procedure  get  next  *) 
done  :=  false  ; 
i  :  =  p-.  r; 
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j  :=  P"-  c; 
casu  temp 
right  : 


left 


up 


down  : 


of 

begin 

j  :=  j+i; 

while  not  done  do  begin 

if  points(.  i-l,j. ). gp  =  k 

then  do_dir(done,up)  (*  go  up  and  done  *) 
else  if (  (subregion(.  k.  ).  list--,  r  =  i) 

and  (subregion(.  k.  ).  list-1,  c  =  j))  then 
begin  (*  return  to  org  point  *) 

new(p); 
p-.  r  :  =  i  ; 

P’-  c  :=  j  ; 
done  :=  true; 

end 

else  if  (points(. i, j+1.  ).  gp  =  k) 

then  j:=  j+1  (*  go  right  *) 

else  if  (points(. i+1, j. ). gp=k  ) 

then  do_dir( done, down)  (*  go  down  and  done  *) 
else  do_dir(done,  left);  (*  go  back  and  done  *) 

end;  (*  while  •') 
end; 
begin 

j  ;=  j-i; 

while  not  done  do  begin 

if  points( . i+1 , j .  ) .  gp  =  k 

then  do_dir(done,down)  (*  go  down  and  done*) 
else  if  ( points (. i, j-1.  ). gp  =  k) 

then  j  :=  j-1  (*  go  left  *) 

else  if  (points( . i-1 , j .  ). gp  =  k) 

then  do_dir(done,up)  (*  go  up  and  done  *) 
else  do_dir(done, right);  (*  go  back  and  done  *) 
end;  (*  while  *) 
end; 
begin 

i  :=  i-1; 

while  not  done  do  begin 

if  points(.  i, j-1.  ).  gp  =  k 

then  do_dir(done, left)  (*  go  left  and  done  *) 
else  if  (points(. i-1, j. )■  gp  =  k) 

then  i  : =  i-1  (*  go  up  *) 

else  if  (points( . i, j  +  i. ).  gp  =  k  ) 

then  do_dir(done,right)(*  go  right  and  done*) 
else  do_dir( done , down);  (*  go  down  and  done  *) 
end;  (*  while  *) 
end; 
begin 

i  :=  i+1; 

while  not  done  do  begin 

if  points(. i, j+1. ). gp  =  k 

then  do_dir( done .right) (*  go  right  and  done*) 
else  if  points(. i+1, j. ). gp  =  k 

then  i  :=  i+1  (*  go  down  *) 

else  if  points( . i , j -1. ) .  gp  =  k 

then  do_dir(done, left)(*  go  left  and  done*) 
else  do_dir(done ,up);  (*  go  up  and  done  *) 
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end;  (*  while  *) 
end; 

end;  (*  case  *) 

end;  (*  sub  procedure  get_next  *) 


begin  (*  procedure  get  boundaries  *) 
for  k  : =  1  to  counter  do 
if  subregion(. k. ).  active  then 
begin 

new(p); 

p-.  r  : =subregion( .  k.  ). upper_left. r  ; 
p-.  c  : =subregion(.  k. ).  upper_left. c  ; 
while  (  points( .  p-.  r-1  ,p-.  c.  ).  gp  =  k)  do 

p-.  r  :=  p-.  r-1;  (*  move  to  the  upper  most  row  of  the  gp  *) 

while  (  points(.  p-.  r,p-.  c-1.  ).  gp  =  k)  do 

p-.  c  :=  p-.  c-1;  (*  move  to  the  left  most  col  of  the  gp  *) 

subregion( . k. ). list  :=  p  ; 
if  points(.  p-.  r,p-.  c+1.  ).  gp  =  k  then 

begin  (*  starting  toward  right  groups  *) 

current  : =  p  ; 
temp  :=  right  ; 
repeat 

get_next(p,temp); 
current-1,  next  :  =  p  ; 
current  : =  p  ; 

until  (subregion!, k. ). list-. r  =  p-. r) 

and  (subregion(. k.  ).  list-. c  =  p-. c); 
if  (temp=down)  and  (points! . p-. r+l,p-. c. ). gp=k)  then 
begin  (*  right  and  down  points  gp  *) 

current  :=  p  ;  (*  hinged  on  start  point  *) 

temp  : =  down  ; 
repeat 

get_next(p , temp) ; 
current-. next  :=  p  ; 
current  : =  p  ; 

until  (subregion(. k. ). list-. r  =  p-. r) 

and  (subregion(.  k.  ).  list-,  c  =  p-.  c); 
p-.  next  :  =  nil; 

end 

else 

p-. next  :=  nil;  (*  right  but  not  down  gp  *) 

end 

else  if  points( . p-. r+1 ,p-. c. ). gp  =  k  then 
begin  (*  starting  toward  down  group  *) 

current  : =  p  ; 
temp  :  =  down  ; 
repeat 

get_next( p , temp) ; 
current-. next  :=  p  ; 
current  : =  p  ; 

until  (subregion(.  k.  ). list-. r  =  p-. r) 

and  (subregion( .  k.  ).  list-. c  =  p-.  c); 
p-.  next  •  =  nil; 
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end 

else 

p-nnext  :=  nil;  (*  one  point  gp  *) 

end; 


the  following  segment  of  code  may  be  used  for  checking  purpose 
for  k  : =  1  to  counter  do 
if  subregion( . k. ). active  then 
begin 

write( ' Boundary  of  group  '  ,k:3,'is  =>'); 
current  :  =  subregion( . k. ) . list; 
while  current  <>  nil  do 
begin 

write( current-1,  r:  2, 1  , '  .current-1,  c:  2,  '=>'  ); 
current  :  =  current-1,  next; 
end; 

writeln; 

end; 


end;  (*  procedure  get_boundaries  ,v) 


*) 


procedure  calc_coeffs; 
var 

k.i, j ,r ,c,gpnum,ptnum  :  integer; 
t,tl  :  point; 

rowchanged.colchanged: boolean; 
begin 

gpnum  :  =  0; 
ptnum  :  =  0; 

for  k  :  =  1  to  counter  do 

if  subregion( . k. ). active  then 
begin 

rowchanged  : =  false; 
colchanged  :=  false; 
for  i  :  =  1  to  3  do 
for  j  : =  1  to  4  do 

subregion( . k. ) . poly_array( . i , j. )  :=  0; 
subregion( .  k.  ).  ee  :=  0; 
subregion( . k. ). ex  :=  0; 
subregion( . k.  ) .  ey  :=  0; 
ptnum  :=  ptnum  +subregion( . k. ) . nurapts  ; 

(*  num  of  linked  pt  *) 

gpnum  :=  gpnum  +1  ;  (*  num  of  active  gp  *) 

t  :=  subregion( . k. ). start; 
while  t. r  <>  0  do 
begin 

tl  :=  t; 

with  subregion( . k. )  do 

getconst(  t.  r.  t. c, poly_array.ee, ex, ey); 
t:  =  points( .  t.  r , t. c.  ). ptnext  ; 
if  (t. r<>0)  and  (t. r<>tl.r)  then 
rowchanged  :=  true; 
if  (t.c<>0)  and  (t.  c<>tl.c)  the-, 
colchanged  :=  true; 

end; 

with  subregion( . k. )  do 
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begin 

if  rowchanged  and  colchanged  then 

gauss(poly_array,a,b,c)  (*  normal  points  *) 

else  if  not( rowchanged)  and  not( colchanged)  then 
begin 

a  :  =  0; 
b  :  =  0; 

c  :=  poly_array( . 3 ,4. ); 

end  (*  one  point  gp  *) 

else  if  not ( rowchanged)  and  colchanged  then 
begin  (*  horizontal  line  group  *) 

a  :  =  0; 

for  j:=  4  downto  2  do 
poly_array(.  3, j.  )  :  = 

poly_array(.  3, j.  )  *  poly_array(.  2,2.  )  + 
poly_array(.  2, j.  )  *  ( -poly_array( . 3 , 2.  ) ); 
poly_array(.  3,4.  )  :  = 

poly_array( .  3 ,4.  )  / 
poly_array( .  3,3.  ); 

poly_array(. 2,4.  )  :  =  (poly_array(.  2,4.  ) 

poly_array( .  3 ,4.  )  * 
poly_array(.  2,3.  ))  / 

poly_array( .  2,2.  ); 
b  :=  poly_array(.  2,4.  ); 
c  :=  poly_array( .  3 ,4.  ); 

end 

else  if  rowchanged  and  not( colchanged)  then 
begin  (*  vertical  line  group  *) 

b  :  —  0 ; 

for  j:=  4  downto  1  do 
poly_array(.  3, j.  )  :  = 

poly_array(. 3, j. )* 
poly_array( . 1 , 1. )  + 
poly_array( .  1 ,  j.  )* 

( -poly_array( . 3,1. )); 
poly_array( . 3 ,4. )  := 

poly_array( . 3 ,4. )  / 
poly_array(.  3,3. ); 

poly„..array( .  1 ,4.  )  :=  (poly_array(.  1,4.  )- 
poly_array(. 3,4. )* 
poly_array( . 1,3. ))/ 
poly_array( . 1,1. ); 
a  :  =  poly_array( . 1 ,4. ); 
c  :  =  poly_array( . 3,4. ); 
end; 

(*  writeln( '  a,b , c ,=' , a: 10 , '  ' ,b: 10,'  ' ,c: 10)  *) 

end;  (*  with  subregion( . k. )  *) 

end;  (*  if  subregion( . k. ). active  *) 

writeln( ' active  Dts=' ,ptnum: 5 , '  *  active  gpnum  =' ,gpnum:  3); 
end;  (*  procedure  calc_coeffs  *) 


procedure  write_planar_data(points:  point_array); 
var 

r,c  :  integer; 
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tempa, tempb, tempc  :  real; 
begin 

for  r: =  2  to  row-1  do 
for  c:  =  2  to  col-1  do 
begin 

tempa  :=  subregion( . points( . r , c. ). gp. ) . a; 
tempb  :=  subregion( . points( . r , c. ). gp. ) . b; 
tempc  :=  subregion(  .  points( . r ,c.  ). gp.  ). c; 
points(. r,c. ). el  :=  trunc  (tempa*r  +  teropb*c  +  tempc); 
end; 

for  c  : =  1  to  col  do 
for  r  :  =  row  downto  1  do 

writeln(doutl ,points( . r,c. ). el); 

end; 

(*  procedure  write_planar_data  *) 


procedure  write_boundary_data( points: point_array); 
var 

k, r ,c , r 1 , r2 ,cl , c2  :  integer; 
tempa , tempb , tempc  :  real; 
current  :  nodepointer; 
begin 

for  r: =  2  to  row-1  do 
for  c: =  2  to  col-1  do 

points(. r,c. ). el  : =  1  ; 
for  k: =  1  to  counter  do 
if  subregion( . k. ). active  then 
begin 

current  :=  subregion( . k. ) . list  ; 

while  current  <>  nil  do 

begin 

r  1:  =currenf.  r; 
cl:  =currenf.  c; 

if  current-. next  <>  nil  then  begin 
r2: =current-. next-. r; 
c2: =current-. next-,  c; 

end 

else  begin 
r2  :=  rl; 
c2  :  =  cl; 
end; 

if  rl  >  r2  then  begin  r:=  r2;  r2:  =  rl;  rl:=  r;  end; 
if  cl  >  c2  then  begin  c:  =  c2;  c2:=  cl;  cl:=  c;  end; 

tempa  :=  subregion( . k. ). a; 
tempb  :=  subregion( . k. ). b; 
tempc  :=  subregion( . k. ). c; 
if  rl=r2  then  for  c  :=  cl  to  c2  do 

points(. rl,c. ). el  :=  trunc  (tempa*rl  +  tempb*c  +  tempc); 
if  cl=c2  then  for  r  :=  rl  to  r2  do 

points( . r , cl. ) . el  :=  trunc  (tempa*r  +  tempb*cl  +  tempc); 
current  : =  current-. next; 
end; 
end; 

for  c  : =  1  to  col  do 
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for  r  : =  row  downto  1  do 

writeln( dout2 , points ( . r,c. ). el); 
end;  (*  write_boundary_data  *) 

(*  procedure  write_boundary_data  *) 
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(*==== — ===== . =====  == 

Main  Program 


begin 
page; 

initialize  ; 

reset(data, 1  data  ell'); 
reset(datam, 'data  magi'); 
reset(datac, ' data  curl'); 
rewrite(doutl , ' data  outl');  (*  data  after  combine  process  *) 
rewrite(dout2, ' tdata  outl');  (*  data  of  boundary  *) 

for  c  : =  1  to  col  do 
for  r  : =  row  downto  1  do 

readln(data,points( . r,c. ). el); 
for  c  : =  1  to  col  do 
for  r  ; =  row  downto  1  do 
begin 

read ln(dat am, tempmag); 
if  tempmag  <  low_bound  then 

points( . r , c. ) . magtype  :=  level 
else  if  tempmag  >  low_bound  then 
points(. r,c. ). magtype  :=  safe 
else  points( . r,c. ). magtype  :=  unsafe; 
points(. r ,c. ). mag  :=  tempmag; 

end; 

for  c  : =  1  to  col  do 
for  r  : =  row  downto  1  do 

readln(datac .point s( . r ,c. ). cur); 

(*printmagcur*)  (*  print  magnitude  and  curvature  *) 

row  link;  (*  link  points  in  row  and  make  subregions  *) 

printsub(points); (*  print  subregions  created  by  row link  *) 

collink;  (*  link  points  in  row  and  column  and  make  subregions  *) 

calc_coeffs;  (*  calculate  plane  expressions  for  the 

subregions  created  *) 
combine(points .subregion); 

(,v  combine  adjacent  regions  if  they 
have  similar  plane  expressions  *) 
calc_coeffs;  (*  calculate  plane  expressions  for  the 

final  subregions  *) 
get_boundaries( subregion); 

(,v  get  boundaries  of  subregions  in  linked  list  *) 
printsub2( points ) ;  (*  print  subregions  created  by  combine  *) 

write_planar_data( points);  (*  write  out  the  elevation  of  the  model  ,v) 
write_boundary_data(points); (*  write  out  the  elevation  of 

the  boundaries  of  the  subregions  *) 


end. 


APPENDIX  C.  PASCAL  CODE  FOR  SMOOTHING  ELEVATION  DATA 
AND  CALCULATING  GRADIENT 


(*Ss300000*) 

program  calc_slope_curv_smooth( input .output); 


(*******************************************lWr******** ***************** 


Author  =  Yee,  Seung  Hee. 

Date  =  3  June  1988. 

Input  =  1.  Elevation  data  acquired  from  Fort.  Hunter  Liggett 
data  base  (40  by  40  square). 

Output  =  1.  File  of  magnitude  of  gradient  of  the  points. 

2.  File  of  curvature  of  the  points  (optional). 

3.  File  of  smoothed  elevation  data  (optional). 
Computer  =  Waterloo  Pascal  on  IBM  370/3033AP, 

Naval  Postgraduate  School. 


Description  = 

This  program  contains  two  major  procedures: 
calc_s lope_curvature  and  smoothing. 

The  procedure  calc_s lope_curvature  calculates  the  slope 

and  curvature  of  each  points ,  the  procedure  smoothing  processes  the 

elevation  data  to  get  smoothed  elevation  data. 

Original  elevation  data  is  represented  in  feet  and  there  are  40*40 
data  cells  in  one  window.  The  distance  between  each  cell  isl2.5 
meters  ,  so  it  needs  to  be  converted  to  get  a  metric  system. 


**• 


•VfVr 


V.  J.  J.  J.  J.  J. 


) 


const 

distance.unit  =  1;  (*  3.28084*12.5  =  41.0105  *) 

row  =  3; 
col  =  3; 

curv_thresh  =  0. 05;  (*  threshold  for  eigenvalues  of  biquadratic  *) 
upbound  =  0.6;  (*  slope  >  60  %  then  unsafe  slope  *) 

lobound  =  0.  1;  (*  slope  <  10  %  then  level  slope  *) 

type 

elrange  =  0. . 10000; 
pointrec  =  record 

el  :  elrange  ; 
slope:  real; 
cur:  char; 
end; 

point_array=array  ( . 1. . row, 1. . col. )  of  pointrec  ; 

var 

points .spoints  :  point_array  ; 

data.outslope.outcur.smoutel.smoutslope.smoutcur  :  text; 

counter , r , c , i ,num:  integer; 

k2 ,k3 ,k4 ,k5 ,k6 .slope, laml , lara2: real; 


procedure  init(var  pt: point_array), 
var 
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r.c.i.j:  integer; 
begin 

for  r  : =  1  to  row  do 
begin 

pt(.  r,  1.  ).  slope  :=  999; 
pt(.  r, col.  ).  slope  :=  999; 
pt(  .  r ,  1.  ).  cur  :  =  1  u' ; 
pt( .  r ,col.  ).  cur  :=  'u'; 
end; 

for  c  : =  1  to  col  do 
begin 

pt( .  1  ,c.  ).  slope  :=  999; 
pt( . row, c. ). slope  :=  999; 
pt(.  l,c.  ).  cur  :  =  'u' ; 
pt( .  row.c.  ). cur  :=  1  u'; 
end; 


end; 

procedure  calc_slope_curvature( var  pt: point_array); 
var 

r , c , i , j :  integer; 

tempi , temp2 , termx .termy .dtermx .dtermy , 
templam  ,termxx,termyy,termxy .term  : real; 
begin 

for  r: =  2  to  row-1  do 
for  c: =  2  to  col-1  do 
begin 


term  :  - 

=  0.0; 

termx  : 

=  0.0; 

termy  : 

=  0.0; 

termxx 

:=  0.0 

termyy 

:  =  0.  0 

termxy 

:  =  0.  0 

for  i  : 

=-l  to 

for  j  : 

=- 1  to 

begin 

term  :=  term  +  pt( .  r+i,c+j.  ). el; 
termx  :=  termx  +  j*pt( . r+i,c+j. ). el; 
termy  :=  termy  +  i*pt( . r+i,c+j. ). el; 
termxx  :=  termxx  +  j*j*pt(. r+i,c+j. ).  el; 
termyy  :=  termyy  +  i*i*pt(. r+i,c+j. ). el; 
termxy  :=  termxy  +  i*j*pt(. r+i,c+j.  ). el; 
end; 

k2  :=  termx/( 6*distance_unit); 

k3  :=  termy/( 6*distance_unit ) ; 

k4  :=  (termxx/2  -  term/3)/distance_unit; 

k5  :=  termxy/(4*distance_unit)  ; 

k6  :=  (termyy/2  -  term/3)/distance  unit; 

laml  :=  ((2*k4+2*k6)  + 

sqrt( abs( sqr( -2*k4-2*k6) -4*( 2*k4*2*k6-sqr(k5) ) ) ) )/2; 
lam 2  :=  ((2*k4+2*k6)  - 

sqrt( abs( sqr( -2*k4-2*k6 ) -4*( 2*k4*2*k6-sqr(k5 ) ) ) ) ) /2; 
if  abs(laml)  <  abs(lam2)  then 
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"  \ 


begin 

templam  :=  laml; 
laml  : =  lam2; 
lam2  :  =  tetnplam; 

end;  (*  eigenvalues  of  bigger  slopenitude  becomes  laml  *) 
slope  :=  sqrt(sqr(k2)  +  sqr(k3)); 

pt(. r,c. ). slope  : =  slope; 

if  ( laml>curv_thresh)  and  ( lam2>curv_thresh)  then 
pt( . r ,c. ). cur  :=  ' d '  (*  depression  *) 

else  if  ( laml<-curv_thresh)  and  ( lam2<-curv_thresh)  then 
pt(. r,c. ). cur  :=  ' h1  (*  hill  or  peak  *) 

else  if  (abs( laml)<=curv_thresh)  and 
(abs( lam2)<=curv_thresh)  then 
pt(. r ,c. ). cur  :=  1 p 1  (*  planar  *) 

else  if  ( laml<-curv_thresh)  and  ( lam2>curv_thresh)  then 
pt(. r,c. ). cur  :=  's'  (*  saddle  *) 

else  if  ( laml<-curv_thresh)  and 

( abs( lam2)<=curv_thresh)  then 
pt(. r,c. ). cur  :=  'r'  (*  ridge  *) 

else  if  ( abs( lam2)<=curv_thresh)  and 
( laml>curv_thresh)  then 
pt( .  r  ,  c.  ) .  cur  :=  V  (*  valley  *) 

else  if  ( lam2<-curv_thresh)  and  ( laml>curv_thresh)  then 
pt(. r,c. ). cur  :  =  'a'  (*  pass  *) 

end; 

end; 


procedure  smoothing( var  points, spoints: point_array); 
var 

r ,c , i , j , temp.dif  :  integer; 
begin 

for  r  :=  2  to  (row-1)  do 
for  c  : =  2  to  (col-1)  do 
begin 

if  (points(. r,c. ). slope  >  lobound) 
then  begin 

temp  :=  0; 

for  i  :=  r-1  to  r+1  do 

temp  :=  points( . i,c. ). el  +  temp; 
for  j  :=  c-1  to  c+1  do 

: emp  :=  points( .  r , j. ). el  +  temp; 
spoints( . r , c. ) . el  :=  trunc( temp/6) ; 

end 

else  spoints(. r,c. ). el  :=  points( . r , c. ) . el  ; 
end; 

for  r  :  =  1  to  row  do 
begin 

spoints( . r , 1. ) . el  :=  points(. r,l. ). el; 
spoints( . r , col. ) . el  :=  points(. r,col. ). el; 
end; 

for  c  : —  1  to  col  do 
begin 

spoints ( . 1 , c. ) . el  : =  point s( . 1 ,c. ). el; 
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spoints(.  row,c.  ).  el  :=  points! . row , c. ).  el; 
end; 

end;  (*  procedure  smoothing  *) 


procedure  check(var  pt: point_array); 
var 

count , countb , countp , counth , countd  : integer; 

count e, count r, countv, count a  : integer; 
r,c,aa,cc  :  integer; 
begin 

count  :=0; 
countb  : =0; 
countp  : =0; 
counth  :=0; 
countd  : =0; 
count e  : =0; 
count r  : =0; 
countv  : =0; 
counta  : =0; 
aa  : =  0; 
cc  :  =  0; 

for  r  : =  1  to  row  do 
for  c  : =  1  to  col  do 
if  pt(. r,c. ). slope  <  lobound  then 
aa  :  =  aa  +  1 

else  if  (pt( . r,c. ). slope  =  999.0)  then 
countb  : =  countb  +  1 

else  if  (pt(. r,c. ). slope  >  upbound)  then 
cc  :  =  cc  +  1 


else 

case 


end 
writeln! 
writeln! 
writeln! 
writeln! 
writeln! 
writeln! 
writeln( 
writeln! 
writeln! 
writeln( 
writeln! 
end; 


r ,c.  ).  cur  of 


1  1 
p 

countp  : =  countp  +1 

'd' 

countd  :=  countd  +1 

'  h' 

counth  :=  counth  +1 

'  e' 

counte  : =  counte  +1 

'  r ' 

countr  :=  countr  +1 

'  v' 

countv  :=  countv  +1 

'  a' 

counta  :=  counta  +1 

level  points  are'faa); 
unsafe  points  are  ,cc); 
for  safe  points: '); 
count  planar  =' , countp); 
count  hill  =' , counth); 
count  depression  =' , countd); 
count  etc  =' ,counte); 
count  saddle  =' , counts); 
count  ridge  ='  countr); 
count  valley  =  .countv); 
count  pass  =' , counta)*) 
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Main  Program 


(* 


*) 


begin 

reset(data, 'data  ell'); 
rewrite(outslope, 'data  ma^l'); 
rewrite(outcur , ' data  curl1); 
rewrite( smoutel , ' smdata  ell'); 
rewr ite( smouts lope , ' smdata  magi' ) ; 
rewrite(smoutcur, ' smdata  curlT); 

for  c  :=  1  to  col  do 
for  r  :  =  row  dowmto  1  do 

readln(data,points(.  r,c. ). el); 
init(points); 
init(spoints); 

calc_s lope_curvature( points ) ; 

(*  check( points)  *) 
for  c  :  =  1  to  col  do 
for  r  : =  row  downto  1  do 
begin 

writeln(outslope,points(. r,c. ). slope); 
writeln(outcur,points(. r,c. ). cur); 
end; 

for  r:=  1  to  10  do  begin 

smoothing(points , spoints); 
calc_s lope_curvature( spoints ) ; 
check( spoints); 
points  : =  spoints; 
end; 

for  c  : =  1  to  col  do 
for  r  : =  row  downto  1  do 
begin 

writeln(smoutel,spoints(. r,c. ). el); 
writeln(smoutslope,spoints(. r,c. ). slope); 
writeln( smoutcur ,spoints( .  r,c. ). cur); 
end; 

end. 
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